9  Dados Geográficos

Conteúdos

Introdução à análise espacial de dados geográficos.

Os formatos de dados JSON e GEOJSON.

Como obter dados através de um API.

9.1 Manipulação de dados

Geometria em GeoPandas

• Geoseries
• Geometry

Operações sobre a coluna GeoMetry

• simplificar é importar
• facilita a visualização

import geopandas as gpd
import folium

# Simplify: https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.simplify.html
# Exemplo Importar os dados da BGRI2021 
# Caminho para o arquivo GeoPackage
gpk = r"data\geo\BGRI2021_1106.gpkg"

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1106 = gpd.read_file(gpk)
gdf1106_2 = gpd.read_file(gpk)
# Simplificar a geografia para uma precisão de 5 metros
# Experimenta - diferentes valores para ver o efeito na geometria
gdf1106_2['geometry'] = gdf1106_2['geometry'].simplify(tolerance=10)

# ------------------
# Mostrar a localização com Folium
# Obter Centroid
centroid = gdf1106.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude  longitude
center_map = [centroid.y, centroid.x]
# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=12, tiles='OpenStreetMap')

# Adicionar Geografia folium map
# folium.GeoJson constructor
folium.GeoJson(gdf1106).add_to(folium_map)

# Mudar a cor
style_function = lambda x: {'fillColor': '#ffffff', 'color': '#000000'}
folium.GeoJson(gdf1106_2, style_function=style_function).add_to(folium_map)

#  Widget para controloar os diferentes layers:
folium.LayerControl().add_to(folium_map)

folium_map
# Visualizar o GeoDataFrame
#gdf1106.plot(column = 'DTMNFR21',
#              legend = False)

Make this Notebook Trusted to load map: File -> Trust Notebook

# explore() forma fácil de ver uma GeoDataFrame

# Mostrar a geografia do GeoDataFrame
gdf1106_2.explore(column = 'DTMNFR21',
              legend = True,
                  edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

Utilizar a área de cada elemento

#print(gdf1106.info())

# Codigo que mostra como calcular a população por KM2

# Verificar o CRS da GDF 
print (gdf1106.crs)
# Area 1º registo
print('BGRI21:', gdf1106.iloc[10].BGRI2021, 'Area:', gdf1106.iloc[10].geometry.area)

# Adicionar nova coluna 
# Caso GDF está noutra CRS será necessario uma correção: gdf1106['geometry'].to_crs(epsg=3857).area
gdf1106['AREA_KM2'] = gdf1106['geometry'].area / 1000000
gdf1106['INDIV_KM2'] = gdf1106['N_INDIVIDUOS'] / gdf1106['AREA_KM2']

# Mostrar resultado:
print(gdf1106[['DTMNFR21', 'N_INDIVIDUOS','AREA_KM2', 'INDIV_KM2']].head(10))

EPSG:3763
BGRI21: 11065601903 Area: 12351.477088540338
  DTMNFR21  N_INDIVIDUOS  AREA_KM2     INDIV_KM2
0   110656         202.0  0.009658  20915.808095
1   110657         197.0  0.008445  23326.786403
2   110658          21.0  0.002582   8131.932073
3   110658          20.0  0.003078   6498.053191
4   110658          22.0  0.002556   8608.476675
5   110607         178.0  0.070032   2541.699854
6   110656         195.0  0.007235  26953.994420
7   110610          12.0  0.002018   5947.325321
8   110657          42.0  0.004894   8581.337906
9   110610          87.0  0.021791   3992.414803

Mostrar resultado como mapa

# Import packages
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

# Definir figura:
f, ax = plt.subplots(1, figsize=(9, 9))

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 3}

# Generate the choropleth and store the axis
# natural_breaks
gdf1106.plot(column=gdf1106.INDIV_KM2, 
              scheme='quantiles', # natural_breaks, quantiles, equal_interval 
              k=9, 
              cmap='PuBu', 
              legend=True,
              edgecolor = 'None', # sem outline
              legend_kwds  = lgnd_kwds,
              ax = ax)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('População por Km2')
plt.show()

Alterar o tipo de Geometria

#print(gdf1106.info())

import geopandas as gpd
import folium

# Criar um novo GeoDataFrame de pontos
gdf1106_points = gdf1106.copy()

# Neste momento a Geometria ainda é de Polygons:
print (gdf1106_points['geometry'].iloc[0].geom_type)

MultiPolygon

converter para pontos

# Calcular o centróide de cada polígono
gdf1106_points['geometry'] = gdf1106['geometry'].centroid

# Calcular ponto dentro poligono
#gdf1106_points['geometry'] = gdf1106['geometry'].representative_point()

print("Geometria original e nova:",gdf1106.iloc[0]['geometry'].geom_type,gdf1106_points.iloc[0]['geometry'].geom_type)

Geometria original e nova: MultiPolygon Point

visualizar resultado

# Fazer Seleção dos Registos para facilitar visualização
gdf_pnt110657 = gdf1106_points[gdf1106_points['DTMNFR21'] == '110655']
gdf_poly110657 = gdf1106[gdf1106['DTMNFR21'] == '110655']

# --------------------------------------
# Mostrar a localização com Folium
# São muitos dados - visualização é lento
centroid = gdf_pnt110657.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude  longitude
center_map = [centroid.y, centroid.x]
# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=15, tiles='OpenStreetMap')

# Adicionar Geografia folium map
# folium.GeoJson constructor
folium.GeoJson(gdf_pnt110657).add_to(folium_map)
folium.GeoJson(gdf_poly110657).add_to(folium_map)

#  Widget para controloar os diferentes layers:
folium.LayerControl().add_to(folium_map)

folium_map

Make this Notebook Trusted to load map: File -> Trust Notebook

criar nova geografia a partir de um dissolve

import geopandas as gpd

# Codigo que mostra um dissolve das subsecções para freguesias
# Utilizar argumento aggfunc (default = first)

# Alternativa 1: DTMNFR21 passa a ser o index - sem reset_index()
# gdf1106_freg = gdf1106.dissolve(by='DTMNFR21', aggfunc='sum')

# Alternativa 2 reset_index para manter a coluna
gdf1106_freg = gdf1106.dissolve(by='DTMNFR21', aggfunc='sum').reset_index()


# Mostrar Resultado da nova gdf
# print(gdf1106_freg.info())

# De seguido será necessária fazer limpeza e correção das colunas
# Apagar Colunas desnecessários
gdf1106_freg = gdf1106_freg.drop(columns=['BGRI2021','DTMNFRSEC21','SECNUM21','SSNUM21','SECSSNUM21','SUBSECCAO','NUTS1','NUTS2','NUTS3'])

# Mudar os valores das colunas nivel superior a DTMNFR21
gdf1106_freg['DTMN21'] = '1106'
gdf1106_freg['DT21'] = '1106'

gdf1106_freg.head()

	DTMNFR21	geometry	OBJECTID	DT21	DTMN21	N_EDIFICIOS_CLASSICOS	N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ	N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS	N_EDIFICIOS_EXCLUSIV_RESID	N_EDIFICIOS_1_OU_2_PISOS	...	N_INDIVIDUOS_H	N_INDIVIDUOS_M	N_INDIVIDUOS_0_14	N_INDIVIDUOS_15_24	N_INDIVIDUOS_25_64	N_INDIVIDUOS_65_OU_MAIS	SHAPE_Length	SHAPE_Area	AREA_KM2	INDIV_KM2
0	110601	POLYGON ((-93010.757 -106447.914, -93010.685 -...	4597315	1106	1106	2716.0	1862.0	848.0	2635.0	1896.0	...	6566.0	7740.0	1796.0	1268.0	7323.0	3919.0	66058.035548	2.876608e+06	2.876608	1.576230e+06
1	110602	POLYGON ((-92096.453 -107111.798, -92120.721 -...	3632570	1106	1106	1493.0	475.0	1013.0	1366.0	578.0	...	6269.0	7581.0	1774.0	1235.0	7500.0	3341.0	66916.297439	5.074775e+06	5.074775	1.298304e+06
2	110607	POLYGON ((-85689.781 -103386.061, -85713.171 -...	2891282	1106	1106	1730.0	1033.0	694.0	1673.0	1046.0	...	5786.0	6397.0	1351.0	1148.0	6542.0	3142.0	53627.590860	2.483785e+06	2.483785	1.065595e+06
3	110608	POLYGON ((-92823.294 -104458.654, -92953.772 -...	4639555	1106	1106	2493.0	1074.0	1407.0	2274.0	1050.0	...	15808.0	19554.0	4101.0	3193.0	17593.0	10475.0	115431.283511	8.024926e+06	8.024926	2.225803e+06
4	110610	POLYGON ((-90527.955 -104449.513, -90599.907 -...	3504714	1106	1106	2255.0	1413.0	828.0	2164.0	1467.0	...	6989.0	7798.0	1702.0	1413.0	8077.0	3595.0	66140.890825	2.774336e+06	2.774336	1.185843e+06

5 rows × 41 columns

# Mostrar o resultado (a geografia do GeoDataFrame)
gdf1106_freg.explore(column = 'DTMNFR21',
              legend = True,
                  edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

9.1.1 Operações entre Datasets

Spatial join

import geopandas as gpd
import folium

# Importar Paragens de autocarro
gpk = r"data\geo\GPK_CARRIS.gpkg"

# Ler os dados do GeoPackage para um GeoDataFrame
gdfCarris = gpd.read_file(gpk,encoding='utf-8')

print(gdfCarris.head())
# Total de 1983 registos
print(gdfCarris.info())


print('BGRI21:', gdfCarris.iloc[10].other_tags)

# ------------------
# Mostrar a localização 
gdfCarris.explore(legend = True,
                  edgecolor = 'black',
                  marker_type = 'marker')

      osm_id                                          name barrier   highway  \
0   25726357                                Alcântara Mar          bus_stop   
1  256063280                                      Limoeiro                     
2  258155049                             Estação Oriente          bus_stop   
3  376036462  Campo Pequeno - Avenida Defensores de Chaves          bus_stop   
4  387760167                        Estação de Sete Rios          bus_stop   

     ref address is_in place man_made  \
0   4002                                
1  13804                                
2  60207                                
3   1904                                
4   1621                                

                                          other_tags  \
0  "bench"=>"no","bus"=>"yes","mapillary"=>"30186...   
1  "bench"=>"yes","bin"=>"yes","lit"=>"no","netwo...   
2  "bench"=>"yes","bin"=>"yes","bus"=>"yes","loca...   
3  "bench"=>"yes","bin"=>"yes","bus"=>"yes","depa...   
4  "bench"=>"yes","bin"=>"yes","bus"=>"yes","depa...   

                         geometry  
0  POINT (-90492.623 -106643.486)  
1  POINT (-86909.942 -105891.249)  
2   POINT (-84059.365 -99553.293)  
3  POINT (-87907.740 -102365.259)  
4  POINT (-89930.824 -102471.872)  
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 1983 entries, 0 to 1982
Data columns (total 11 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   osm_id      1983 non-null   object  
 1   name        1983 non-null   object  
 2   barrier     1983 non-null   object  
 3   highway     1983 non-null   object  
 4   ref         1983 non-null   object  
 5   address     1983 non-null   object  
 6   is_in       1983 non-null   object  
 7   place       1983 non-null   object  
 8   man_made    1983 non-null   object  
 9   other_tags  1983 non-null   object  
 10  geometry    1983 non-null   geometry
dtypes: geometry(1), object(10)
memory usage: 170.5+ KB
None
BGRI21: "bench"=>"yes","bus"=>"yes","network"=>"Carris Metropolitana","network:wikidata"=>"Q111611112","official_name"=>"Pç. Areeiro","public_transport"=>"platform","shelter"=>"yes"

Make this Notebook Trusted to load map: File -> Trust Notebook

Spatial join de polygons para points

# Realizar o spatial join
gdf_join = gpd.sjoin(gdf1106_freg, gdfCarris, how="inner", predicate="contains")

# Limitar Colunas:
# gdf_join = gpd.sjoin(gdf1106_freg[['DTMNFR21', 'geometry']], gdfCarris, how="left", predicate="contains")

print (f"Tipo de Geometria: {gdf_join['geometry'].iloc[0].geom_type}")
print (f"Nº de Registos Input: {len(gdf1106_freg)}")
print (f"Nº de Registos Output: {len(gdf_join)}", "\n")


#gdf_join.info()

Tipo de Geometria: Polygon
Nº de Registos Input: 24
Nº de Registos Output: 1983 

Join de points para polygons

por exemplo, obter a freguesia para cada paragem de autocarro

import geopandas as gpd

# Perform spatial join
gdf_join = gpd.sjoin(gdfCarris, gdf1106_freg[['DTMNFR21', 'geometry']], how='left', predicate='within')

# Mostrar Resultado
print (f"Tipo de Geometria: {gdf_join['geometry'].iloc[0].geom_type}")
print (f"Nº de Registos Input: {len(gdfCarris)}")
print (f"Nº de Registos Output: {len(gdf_join)}", "\n")

gdf_join[['osm_id','DTMNFR21']].head()

Tipo de Geometria: Point
Nº de Registos Input: 1983
Nº de Registos Output: 1983 

	osm_id	DTMNFR21
0	25726357	110660
1	256063280	110665
2	258155049	110662
3	376036462	110657
4	387760167	110639

obter a contagem das paragens de freguesia

#gdf1106_freg.drop(columns=['n_paragens_x','n_paragens_y'], inplace=True)

# Realiza o spatial join
# '''
# 1. Fazer sjoin: resultado um GDF com o memso numero de registos que os pontos
# 2. Adicionar uma nova coluna n_paragens com total de registos existentes por Freguesia
# 3. Obter Dataframe com numero de Valores unicos 
# 4. Fazer merge do novo valor obtido com Geodataframe original
# 5. Apagar o objecto do join
# '''
# Realizar o spatial join
# Resultado terá o mesmo nº de registos que gdfCarris
gdf_join = gpd.sjoin(gdf1106_freg, gdfCarris, how="inner", predicate="contains")

print (f"Nº de Registos Output Join: {len(gdf_join)}")


# Contar o número de pontos em cada polígono (novo atributo n_paragens)
# Informação está duplicada para cada Freguesia
gdf_join["n_paragens"] = gdf_join.groupby("DTMNFR21")["geometry"].transform("size")

# # Selecionar o primeiro valor de 'n_paragens' dentro de cada grupo 'DTMNFR21'
unique_values = gdf_join.groupby('DTMNFR21')['n_paragens'].first().reset_index()

print(unique_values.info())

# Exibir o DataFrame resultante
#print(unique_values.head())

# Apagar coluna - caso repetir o codigo (o merge) sem recriar gdf1106
# gdf1106_freg.drop(columns=['n_paragens','n_paragens_x','n_paragens_y'], inplace=True)
gdf1106_freg = gdf1106_freg.merge(unique_values, on='DTMNFR21', how='left')

del gdf_join


# Exibir a GeoDataFrame resultante com o atributo 'n_paragens' adicionado
# print(gdf1106_freg.info())
print(gdf1106_freg[['DTMNFR21', 'n_paragens']].head(10))

Nº de Registos Output Join: 1983
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 24 entries, 0 to 23
Data columns (total 2 columns):
 #   Column      Non-Null Count  Dtype 
---  ------      --------------  ----- 
 0   DTMNFR21    24 non-null     object
 1   n_paragens  24 non-null     int64 
dtypes: int64(1), object(1)
memory usage: 516.0+ bytes
None
  DTMNFR21  n_paragens
0   110601          71
1   110602          72
2   110607          53
3   110608         132
4   110610          75
5   110611          73
6   110618         137
7   110621         138
8   110633         134
9   110639          85

# Desenhar mapa com resultado:
import geopandas as gpd
import matplotlib.pyplot as plt

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 2}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf1106_freg.plot(column=gdf1106_freg.n_paragens, 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=5, 
                      cmap='YlGn', 
                      legend=True,
                      edgecolor = 'None', 
                      legend_kwds  = lgnd_kwds)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Nº de paragens Carris')
plt.show()

outra possibilidade - fazer a seleção de nº de pontos para freguesia específica

# Fazer seleção de nº de pontos para freguesia especifica
from geopandas.tools import sjoin

# Selecionar o polígono específico da freguesia 110655 (Areeiro)
poligono_especifico = gdf1106_freg[gdf1106_freg['DTMNFR21'] == '110655']

# Fazer Join com gdfCArris para obter os pontos
joined = sjoin(gdfCarris, poligono_especifico, how='inner', predicate='within')

# Contar quantos pontos estão dentro do polígono específico
quantidade_pontos = len(joined)

# Mostrar o resultado
print(f"Quantidade de paragens de autocarro dentro a freguesia 110655: {quantidade_pontos}")

Quantidade de paragens de autocarro dentro a freguesia 110655: 49

efectuar um overlay

exemplo de cortar um círculo (buffer) em volta e um ponto

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import Point

# Criar variáveis para a figura
f, ax = plt.subplots(1, figsize=(9, 9))

# Criar um objeto Point
# Atenção primeiro o valor x (longitude) e depois o y (latitutde) - Google Maps devolve Latitude, Longitude
ponto = Point(-9.184111016,38.768216306)
# 38.76821630632057, -9.184111016081756
# Cria um GeoDataFrame do ponto com CRS WGS84
gdf_ponto = gpd.GeoDataFrame([1], geometry=[ponto], crs='EPSG:4326')

# Mudar a projeção do pontos para a projeção da gdf1106:
gdf_ponto = gdf_ponto.to_crs(gdf1106.crs)

print (f"Tipo de Geometria: {gdf_ponto['geometry'].iloc[0].geom_type}")

# Cria um buffer de 500 metros em volta do ponto
gdf_ponto['geometry'] = gdf_ponto.geometry.buffer(1500)
print (f"Tipo de Geometria: {gdf_ponto['geometry'].iloc[0].geom_type}")

# Realiza a interseção entre o buffer e os polígonos
# Opções: intersection’, ‘union’, ‘identity’, ‘symmetric_difference’ or ‘difference’ 
intersecao = gpd.overlay(gdf1106, gdf_ponto, how='symmetric_difference')

# Visualizar o Resultado
intersecao.plot(column = 'DTMNFR21',
              legend = False,
               ax = ax)

# Add basemap do contextily
ctx.add_basemap(
    ax,
    crs=intersecao.crs,
    source=ctx.providers.CartoDB.VoyagerNoLabels,
)


ax.set_axis_off()

plt.show()

Tipo de Geometria: Point
Tipo de Geometria: Polygon

9.1.2 Exportar dados GeoDataframe

principais formatos de output

• Shapefile: Sem necessidade de especificar driver
• GeoJSON: driver=‘GeoJSON’
• GeoPackage: driver=‘GPKG’

gdf1106_freg.to_file(r'data\geo\c2021_fr1106.gpkg', layer='FR1106', driver="GPKG")

listar outputs possíveis

import fiona
fiona.supported_drivers

{'DXF': 'rw',
 'CSV': 'raw',
 'OpenFileGDB': 'raw',
 'ESRIJSON': 'r',
 'ESRI Shapefile': 'raw',
 'FlatGeobuf': 'raw',
 'GeoJSON': 'raw',
 'GeoJSONSeq': 'raw',
 'GPKG': 'raw',
 'GML': 'rw',
 'OGR_GMT': 'rw',
 'GPX': 'rw',
 'Idrisi': 'r',
 'MapInfo File': 'raw',
 'DGN': 'raw',
 'PCIDSK': 'raw',
 'OGR_PDS': 'r',
 'S57': 'r',
 'SQLite': 'raw',
 'TopoJSON': 'r'}

exemplo para obter informação a partir de uma seleção dentro dum buffer

importar dados

import geopandas as gpd
import matplotlib.pyplot as plt

# Simplify: https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.simplify.html
# Exemplo Importar os dados da BGRI2021 
# Caminho para o arquivo GeoPackage
gpk = r'data\geo\BGRI2021_1106.gpkg'

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1106 = gpd.read_file(gpk,encoding='utf-8')

mostrar mapa como follium

# Function para mostrar folium
def mostrarFolium(gdfExtent,gdfDesenhar):
    centroid = gdfExtent.to_crs(epsg=4326).unary_union.centroid

    # Criar Listagem com localização de latitude  longitude
    center_map = [centroid.y, centroid.x]
    # Criar Mapa e mostrar
    folium_map = folium.Map(location=center_map, zoom_start=15, tiles='OpenStreetMap')

    # Adicionar Geografia folium map
    # folium.GeoJson constructor
    for lay in gdfDesenhar:
        folium.GeoJson(lay).add_to(folium_map)
    
    folium_map

executar o buffer

from shapely.geometry import Point

# Cria um objeto Point (podem escolher um ponto no google maps)
# Atenção primeiro o valor x (longitude) e depois o y (latitutde)
ponto = Point(-9.137616,38.738561)

# Cria um GeoDataFrame do ponto com CRS WGS84
gdf_ponto = gpd.GeoDataFrame([1], geometry=[ponto], crs='EPSG:4326')

# Manteer ponto Original
gdf_ponto_original = gdf_ponto.copy()

# Mudar a projeção:
gdf_ponto = gdf_ponto.to_crs(gdf1106.crs)

# Fazer o Buffer
gdf_ponto['geometry'] = gdf_ponto.geometry.buffer(500)

# Executar o Spatial join entre o ponto e as subsecções da BGRI2021
join = gpd.sjoin(gdf1106, gdf_ponto, how="inner", predicate='intersects')

# Somar os valores do atributo N_INDIVIDUOS
soma = join['N_INDIVIDUOS'].sum()

print(f"Soma de N_INDIVIDUOS nos polígonos selecionados: {soma} em volta do ponto ({ponto.y},{ponto.x})")


#-----------------------------------------------------------------
# Mostrar 
mostrarFolium(gdf_ponto,[gdf_ponto,gdf_ponto_original])

Soma de N_INDIVIDUOS nos polígonos selecionados: 13943.0 em volta do ponto (38.738561,-9.137616)

9.2 Desafio

• Definir um ponto, valores latitude e longitude
• Criar um Buffer em volta desse ponto
• Fazer o cálculo da população nas subsecções envolventes

Podem fazer download de outros GPK no site do INE: https://mapas.ine.pt/download/index2021.phtml

# Exemplo Importar os dados da BGRI2021 para o municipio de Valongo
# Caminho para o arquivo GeoPackage
gpk = r'data\geo\BGRI2021_1315.gpkg'

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1315 = gpd.read_file(gpk,encoding='utf-8')

from shapely.geometry import Point

# Cria um objeto Point (podem escolher um ponto no google maps)
# Atenção primeiro o valor x (longitude) e depois o y (latitutde)
# ponto do google maps Sta Rita (Ermesinde)
# 41.20722032311807, -8.541960984473912
ponto = Point(-8.541960984473912,41.20722032311807)

# Cria um GeoDataFrame do ponto com CRS WGS84
gdf_ponto = gpd.GeoDataFrame([1], geometry=[ponto], crs='EPSG:4326')

# Manteer ponto Original
gdf_ponto_original = gdf_ponto.copy()

# Mudar a projeção:
gdf_ponto = gdf_ponto.to_crs(gdf1315.crs)

# Fazer o Buffer
gdf_ponto['geometry'] = gdf_ponto.geometry.buffer(500)

# Executar o Spatial join entre o ponto e as subsecções da BGRI2021
join = gpd.sjoin(gdf1315, gdf_ponto, how="inner", predicate='intersects')

# Somar os valores do atributo N_INDIVIDUOS
soma = join['N_INDIVIDUOS'].sum()

print(f"Soma de N_INDIVIDUOS nos polígonos selecionados: {soma} em volta do ponto ({ponto.y},{ponto.x})")

mostrarFolium(gdf_ponto,[gdf_ponto,gdf_ponto_original])

Soma de N_INDIVIDUOS nos polígonos selecionados: 4981.0 em volta do ponto (41.20722032311807,-8.541960984473912)

9.3 GeoEstatistica

• Análise de Padroes espaciais

  • autocorrelação espacial
  • Média do vizinho mais próximo

• Mapeamento de clusters

  • Análise de clusters e outliers
  • Análise de Hot Spots
  • Clusterização multivariada
  • Indices compostos

• Modelação de relações espaciais

  • Regressão linear
  • Regressão geograficamente ponderada
  • Regressão multiescalar geograficamente opnderada
  • Minimos quadrados
  • Relações locais bivariadas

9.3.1 Autocorrelação Espacial

Clacular medidas

- Moran Global
- Moran Local
- Getis and Ord's local statistics

import matplotlib.pyplot as plt  # Graphics
from matplotlib import colors
import seaborn as sns  # Graphics
import geopandas as gpd # Spatial data manipulation
import pandas as pd  # Tabular data manipulation
# Para Evitar Aviso Point Patterns
from shapely.geometry import Point

# Bibliotecas pysal
import pysal.lib # importação geral
from pysal.explore import esda  # Exploratory Spatial analytics
from pysal.lib import weights  # Spatial weights
import contextily  # Background tiles

# Bibliotecas última parte notebook exemplo
#import rioxarray  # Surface data manipulation
#import xarray  # Surface data manipulation

importar dados do Geopackage com variáveis do C2021 da BRGI2021

import geopandas as gpd
import matplotlib.pyplot as plt

# Caminho para o arquivo GeoPackage
gpk = r'data\geo\BGRI2021_1106.gpkg'

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1106 = gpd.read_file(gpk)

# Simplificar a geografia para uma precisão de 5 metros
gdf1106['geometry'] = gdf1106['geometry'].simplify(tolerance=5)

# Visualizar o GeoDataFrame
gdf1106.plot(column = 'DTMNFR21',
              legend = False)

<Axes: >

adicionar alguns atributos que permitem fazer a análise

# Calcular Novo Atributo Racio de População 65+ anos
gdf1106['IND65'] = gdf1106.N_INDIVIDUOS_65_OU_MAIS/gdf1106.N_INDIVIDUOS

# Calcular Outros Atributos de interesse para Analisar: , N_EDIFICIOS_3_OU_MAIS_PISOS, N_INDIVIDUOS_H, N_INDIVIDUOS_M
gdf1106['IND14'] = gdf1106.N_INDIVIDUOS_0_14/gdf1106.N_INDIVIDUOS
gdf1106['IND_H'] = gdf1106.N_INDIVIDUOS_H/gdf1106.N_INDIVIDUOS
gdf1106['IND_M'] = gdf1106.N_INDIVIDUOS_M/gdf1106.N_INDIVIDUOS
gdf1106['EDIF_3PISOS'] = gdf1106.N_EDIFICIOS_3_OU_MAIS_PISOS/gdf1106.N_EDIFICIOS_CLASSICOS


# Mostrar Dados
print(gdf1106[['BGRI2021', 'DTMNFR21', 'N_INDIVIDUOS_65_OU_MAIS', 'N_INDIVIDUOS', 'IND65', 'IND14', 'IND_H', 'IND_M','EDIF_3PISOS']].head(10))


# Manter apenas as colunas de interesse: (não é necessário - simplifica o GeoDataFrame)
manter_colunas = ['geometry','BGRI2021', 'DTMNFR21', 'N_INDIVIDUOS_65_OU_MAIS', 'N_INDIVIDUOS', 'IND65','IND14', 'IND_H', 'IND_M','N_EDIFICIOS_CLASSICOS','EDIF_3PISOS']
gdf1106 = gdf1106.loc[:, manter_colunas]

#print(gdf1106.info())

      BGRI2021 DTMNFR21  N_INDIVIDUOS_65_OU_MAIS  N_INDIVIDUOS     IND65  \
0  11065602301   110656                     30.0         202.0  0.148515   
1  11065700203   110657                     62.0         197.0  0.314721   
2  11065801011   110658                      3.0          21.0  0.142857   
3  11065801012   110658                      7.0          20.0  0.350000   
4  11065801013   110658                      4.0          22.0  0.181818   
5  11060701205   110607                     28.0         178.0  0.157303   
6  11065602601   110656                     39.0         195.0  0.200000   
7  11061000601   110610                      1.0          12.0  0.083333   
8  11065700110   110657                     12.0          42.0  0.285714   
9  11061001305   110610                     19.0          87.0  0.218391   

      IND14     IND_H     IND_M  EDIF_3PISOS  
0  0.069307  0.524752  0.475248     0.933333  
1  0.131980  0.380711  0.619289     1.000000  
2  0.190476  0.571429  0.428571     0.250000  
3  0.200000  0.400000  0.600000     0.000000  
4  0.318182  0.409091  0.590909     0.000000  
5  0.095506  0.516854  0.483146     0.162162  
6  0.117949  0.420513  0.579487     0.785714  
7  0.000000  0.416667  0.583333     0.000000  
8  0.047619  0.404762  0.595238     1.000000  
9  0.160920  0.505747  0.494253     0.714286  

tartar NaN

# Contar o número total de NaNs no DataFrame
total_nans = gdf1106.isna().sum().sum()
print('Número total de registros com NaN:', total_nans)

# Contar o número de NaNs em cada coluna
nans_por_coluna = gdf1106.isna().sum()
print('Número de registros com NaN por coluna:\n', nans_por_coluna)

# Corrigir NaN
# Existem 2 possibilidades
# 1. Deixar fora
#gdf1106 = gdf1106.dropna()
# 2. Substituir por outros valores
gdf1106 = gdf1106.fillna(0)

Número total de registros com NaN: 781
Número de registros com NaN por coluna:
 geometry                     0
BGRI2021                     0
DTMNFR21                     0
N_INDIVIDUOS_65_OU_MAIS      0
N_INDIVIDUOS                 0
IND65                      157
IND14                      157
IND_H                      157
IND_M                      157
N_EDIFICIOS_CLASSICOS        0
EDIF_3PISOS                153
dtype: int64

visualização inicial

# Set up figure and a single axis
f, ax = plt.subplots(1, figsize=(9, 9))
# Build choropleth
gdf1106.plot(
    column="IND65",
    cmap="viridis",
    scheme="quantiles",
    k=5,
    edgecolor="white",
    linewidth=0.0,
    alpha=0.75,
    legend=True,
    legend_kwds=dict(loc=2),
    ax=ax,
)
# Add basemap
contextily.add_basemap(
    ax,
    crs=gdf1106.crs,
    source=contextily.providers.CartoDB.VoyagerNoLabels,
)
# Remove axes
ax.set_axis_off()

plt.show()

9.3.2 Matriz de vizinhança e Moran Global

Queen

f we wanted them to be considered as neighbours, we can switch to the more inclusive notion of Queen contiguity, which requires the pair of polygons to only share one or more vertices. We can create the neighbor relations for this same configuration as follows:

KNN

The first type of distance based weights defines the neighbor set of a particular observation as containing its nearest observations, where the user specifies the value of . To illustrate this for the San Diego tracts, we take . This still leaves the issue of how to measure the distance between these polygon objects, however. To do so we develop a representative point for each of the polygons using the centroid.

from pysal.model import spreg

# Input gdf1106 (BGRI de Lisboa)

# Calcular a matriz de pesos espaciais
# Metodo Original: 0.13069266093679838 e Valor-p: 0.001
w = pysal.lib.weights.Queen.from_dataframe(gdf1106, use_index=True)
# Metodo knn (I de Moran: 0.11720923169446687; Valor-p: 0.001)
#w = pysal.lib.weights.KNN.from_dataframe(gdf1106, k=5)

# Lidar com ilhas
# Normalizar a matriz de pesos
# Row-standardization 
w.transform = "r"

# Corrigir para ilhas
#w.set_transform('r')
islands = w.islands
if islands:
    for island in islands:
        w.neighbors[island] = [island]
        w.weights[island] = [1]  

# Selecionar a coluna com a percentagem de população superior a 65 anos
data = gdf1106['IND65']

# Calcular a estatística de Moran
moran = esda.Moran(data, w)

# Imprimir o valor I de Moran e o valor-p
print('I de Moran:', moran.I)
print('Valor-p:', moran.p_sim)

# Plotar o gráfico de dispersão de Moran (desatividado - vai ser efetuado em tarefas na próxima seção)
#plot_moran(moran, zstandard=True, fill=True, figsize=(10,4))

I de Moran: 0.15343960749795466
Valor-p: 0.001

Interpretação do valor Moran Global

\(H_0: I = 0\) (padran aleatorio) vs \(H_1 \ne 0\) (padrão espacial com clusters)

9.3.3 Motivating local spatial autocorrelation (Moran Plot)

Nesta Parte é efetuado o seguinte:

• Calcular o Spatial Lag
• Calcular as versões centradas
• Criar um Scatterplot
• Visualizar a distribuição dos valores em relação ao médio e aos valores nos polígonos vizinhos

Moran Plot

The Moran Plot is a way of visualizing a spatial dataset to explore the nature and strength of spatial autocorrelation. It is essentially a traditional scatter plot in which the variable of interest is displayed against its spatial lag. In order to be able to interpret values as above or below the mean, the variable of interest is usually standardized by subtracting its mean:

Calcular o Spatial log

# Calcular o Spatial Lag
# 
gdf1106["w_IND65"] = weights.lag_spatial(w, gdf1106['IND65'])

# And their respective centered versions, where we subtract the average off of every value
gdf1106["IND65_std"] = gdf1106["IND65"] - gdf1106["IND65"].mean()
gdf1106["w_IND65_std"] = weights.lag_spatial(w, gdf1106['IND65_std'])

# Visualizar Valores em Cima e Baixo do Médio
# Set up the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values

sns.regplot(
    x="IND65_std", y="w_IND65_std", data=gdf1106, ci=None
)
plt.show()

Figura com a relação das vizinhanças (Moran Plot)

Mesma figura com indicação dos 4 quadrantos

• LH: Valores na subsecção em baixo do médio, valores circundantes em cima do médio
• HH: Valores na subsecção em cima do médio, valores circundantes em cima do médio
• LL: Valores na subsecção em baixo do médio, valores circundantes em baixo do médio
• HL: Valores na subsecção em cima do médio, valores circundantes em baixo do médio

# Criar os Quadrantos
# Set up the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values
sns.regplot(
    x="IND65_std", y="w_IND65_std", data=gdf1106, ci=None
)

# Esta Parte demora muito tempo

# Add vertical and horizontal lines (definição valor onde adicionar)
plt.axvline(0, c="k", alpha=0.5)
plt.axhline(0, c="k", alpha=0.5)

# Adicionar Text para Cada Quadrant - coordinados Text tendo em conta a distribuição dos dados
plt.text(0.6, 0.20, "HH", fontsize=25, c="r")
plt.text(0.6, -0.15, "HL", fontsize=25, c="r")
plt.text(-0.2, 0.20, "LH", fontsize=25, c="r")
plt.text(-0.15, -0.15, "LL", fontsize=25, c="r")

# Displaby
plt.show()

9.3.4 Local Moran’s I

Calcular LISA (Local Indicators of Spatial Association) e Visualização Inicial

# https://pysal.org/esda/generated/esda.Moran_Local.html
# lisa = moran_loc
# Diferença com Notebook da Formação - dá erro porque os valores são NaN
from splot.esda import lisa_cluster
# data = 
lisa = esda.moran.Moran_Local(gdf1106['IND65'], w)

# Plotar o mapa de clusters de Moran
lisa_cluster(lisa, gdf1106, figsize=(9,9))

(<Figure size 864x864 with 1 Axes>, <Axes: >)

Kernel Estimate Plotting

import numpy as np

# Valores LISa primeiros 10 registos
print(lisa.Is[:10])

# alguns indacores dos valores
print(f'''Minimum: {np.min(lisa.Is)}
Maximum: {np.max(lisa.Is)}
STD: {np.std(lisa.Is)}''')

[ 0.15105448 -0.01156592  0.45131911  0.98016347  0.13006217 -0.22294648
  0.04772319 -0.98709029  0.0033492   0.01171758]
Minimum: -6.358580055279035
Maximum: 5.47984335779404
STD: 0.6495277258099246

import warnings

# Nao mostrar aviso FutereWarning (não aconselhável)
warnings.filterwarnings("ignore", category=FutureWarning)

# Draw KDE line
ax = sns.kdeplot(lisa.Is)

# Add one small bar (rug) for each observation
# along horizontal axis
sns.rugplot(lisa.Is, ax=ax)

plt.show()

Visualizações diferentes, mapas com medidas LISA

Significado dos 4 mapas: 1. Valor LISA cada polígono (valor lisa.Is) 2. Valor do quadrante cada área (esdaplot.lisa_cluster, p = 1) 3. Indicação da significância estatística (lisa.p_sim < 0.05) 4. Combinação dos anterior 2 (esdaplot.lisa_cluster, p = 0.05)

from splot import esda as esdaplot


# Set up figure and axes
f, axs = plt.subplots(nrows=2, ncols=2, figsize=(12, 12))
# Make the axes accessible with single indexing
axs = axs.flatten()

# Subplot 1 #
# Choropleth of local statistics
# Grab first axis in the figure
ax = axs[0]
# Assign new column with local statistics on-the-fly
gdf1106.assign(
    Is=lisa.Is
    # Plot choropleth of local statistics
).plot(
    column="Is",
    cmap="plasma",
    scheme="quantiles",
    k=5,
    edgecolor="white",
    linewidth=0.1,
    alpha=0.75,
    legend=True,
    ax=ax,
)

# Subplot 2 #
# Quadrant categories
# Grab second axis of local statistics
ax = axs[1]
# Plot Quadrant colors (note to ensure all polygons are assigned a
# quadrant, we "trick" the function by setting significance level to
# 1 so all observations are treated as "significant" and thus assigned
# a quadrant color
esdaplot.lisa_cluster(lisa, gdf1106, p=1, ax=ax)

# Subplot 3 #
# Significance map
# Grab third axis of local statistics
ax = axs[2]
#
# Find out significant observations
labels = pd.Series(
    1 * (lisa.p_sim < 0.05),  # Assign 1 if significant, 0 otherwise
    index=gdf1106.index  # Use the index in the original data
    # Recode 1 to "Significant and 0 to "Non-significant"
).map({1: "Significant", 0: "Non-Significant"})
# Assign labels to `gdf1106` on the fly
gdf1106.assign(
    cl=labels
    # Plot choropleth of (non-)significant areas
).plot(
    column="cl",
    categorical=True,
    k=2,
    cmap="Paired",
    linewidth=0.1,
    edgecolor="white",
    legend=True,
    ax=ax,
)


# Subplot 4 #
# Cluster map
# Grab second axis of local statistics
ax = axs[3]
# Plot Quadrant colors In this case, we use a 5% significance
# level to select polygons as part of statistically significant
# clusters
esdaplot.lisa_cluster(lisa, gdf1106, p=0.05, ax=ax)

# Figure styling #
# Set title to each subplot
for i, ax in enumerate(axs.flatten()):
    ax.set_axis_off()
    ax.set_title(
        [
            "Local Statistics",
            "Scatterplot Quadrant",
            "Statistical Significance",
            "Moran Cluster Map",
        ][i],
        y=0,
    )
# Tight layout to minimize in-between white space
f.tight_layout()

# Display the figure
plt.show()

Diferentes contagens (Atributo “q”)

• Atributo ‘q’ do Moran Local mostra os valores nos 4 quadrantos

# Criar pd.series do resultado
counts = pd.Series(lisa.q).value_counts() # Original: pd.value_counts(lisa.q)
counts

# Visualizar os valores do primeiros 10 registos 

lisa.q[:10]

# Visualizar o numero de áreas com valores estatitisticamente significante

(lisa.p_sim < 0.05).sum() * 100 / len(lisa.p_sim)

16.087880935506732

9.3.5 Getis and Ord’s local statistics

Outras medidas para representar a autocorrelação espacial (Getis and Ord’s)

• the Gⁱ: statistic, which omits the value at a site in its local summary
• the G^i*: statistic, which includes the site’s own value in the local summary

calcular indicadores

# Gi
go_i = esda.getisord.G_Local(gdf1106["IND65"], w) # , star=None
# Gi*
go_i_star = esda.getisord.G_Local(gdf1106["IND65"], w, star=True)

C:\Users\bruno.lima\AppData\Roaming\Python\Python312\site-packages\esda\getisord.py:615: UserWarning: Gi* requested, but (a) weights are already row-standardized, (b) no weights are on the diagonal, and (c) no default value supplied to star. Assuming that the self-weight is equivalent to the maximum weight in the row. To use a different default (like, .5), set `star=.5`, or use libpysal.weights.fill_diagonal() to set the diagonal values of your weights matrix and use `star=None` in Gi_Local.
  warnings.warn(

visualizar o resultado

def g_map(g, db, ax, title):
    """
    Create a cluster map
    ...

    Arguments
    ---------
    g      : G_Local
             Object from the computation of the G statistic
    db     : GeoDataFrame
             Table aligned with values in `g` and containing
             the geometries to plot
    ax     : AxesSubplot
             `matplotlib` axis to draw the map on

    Returns
    -------
    ax     : AxesSubplot
             Axis with the map drawn
    """
    ec = "0.8"

    # Break observations into significant or not
    sig = g.p_sim < 0.05

    # Plot non-significant clusters
    ns = db.loc[sig == False, "geometry"]
    ns.plot(ax=ax, color="lightgrey", edgecolor=ec, linewidth=0.1)
    # Plot HH clusters
    hh = db.loc[(g.Zs > 0) & (sig == True), "geometry"]
    hh.plot(ax=ax, color="red", edgecolor=ec, linewidth=0.1)
    # Plot LL clusters
    ll = db.loc[(g.Zs < 0) & (sig == True), "geometry"]
    ll.plot(ax=ax, color="blue", edgecolor=ec, linewidth=0.1)
    # Style and draw
    contextily.add_basemap(
        ax,
        crs=db.crs,
        source=contextily.providers.Esri.WorldTerrain,
    )
    # Flag to add a star to the title if it's G_i*
    st = ""
    if g.star:
        st = "*"
    # Add title
    ax.set_title(f"G{st} statistic {title}", size=15)
    # Remove axis for aesthetics
    ax.set_axis_off()
    return ax

# Set up figure and axes
f, axs = plt.subplots(1, 2, figsize=(12, 6))
# Loop over the two statistics
for g, ax in zip([go_i, go_i_star], axs.flatten()):
    # Generate the statistic's map
    ax = g_map(g, gdf1106, ax, 'Rácio Popuação 65+ Anos')
# Tight layout to minimise blank spaces
f.tight_layout()
# Render
plt.show()

9.3.6 Exercício

**Repetir o código com outras variáveis ecom município de Porto ou outro município

• Utilizar o codigo Anterior para Calcular para outras variaveis
  • Lista variáveis: ‘IND65’,‘IND14’, ‘IND_H’, ‘IND_M’,‘EDIF_3PISOS’
• Repete a mesma análise para outro municipio
  • Faz download do gpk de um municipio no link: https://mapas.ine.pt/download/index2021.phtml

# Caminho para o arquivo GeoPackage do municipio de Valongo
gpk = r"data\geo\BGRI2021_1315.gpkg"

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1315 = gpd.read_file(gpk)

# Calcular Novo Atributo Racio de População 65+ anos
gdf1315['IND65'] = gdf1315.N_INDIVIDUOS_65_OU_MAIS/gdf1315.N_INDIVIDUOS

# Calcular Outros Atributos de interesse para Analisar: , N_EDIFICIOS_3_OU_MAIS_PISOS, N_INDIVIDUOS_H, N_INDIVIDUOS_M
gdf1315['IND14'] = gdf1315.N_INDIVIDUOS_0_14/gdf1315.N_INDIVIDUOS
gdf1315['IND_H'] = gdf1315.N_INDIVIDUOS_H/gdf1315.N_INDIVIDUOS
gdf1315['IND_M'] = gdf1315.N_INDIVIDUOS_M/gdf1315.N_INDIVIDUOS
gdf1315['EDIF_3PISOS'] = gdf1315.N_EDIFICIOS_3_OU_MAIS_PISOS/gdf1315.N_EDIFICIOS_CLASSICOS


# Mostrar Dados
print(gdf1315[['BGRI2021', 'DTMNFR21', 'N_INDIVIDUOS_65_OU_MAIS', 'N_INDIVIDUOS', 'IND65', 'IND14', 'IND_H', 'IND_M','EDIF_3PISOS']].head(10))


from pysal.model import spreg

# Input gdf1106 (BGRI de Lisboa)

# Calcular a matriz de pesos espaciais
# Metodo Original: 0.13069266093679838 e Valor-p: 0.001
w = pysal.lib.weights.Queen.from_dataframe(gdf1315, use_index=True)
# Metodo knn (I de Moran: 0.11720923169446687; Valor-p: 0.001)
#w = pysal.lib.weights.KNN.from_dataframe(gdf1106, k=5)

# Lidar com ilhas
# Normalizar a matriz de pesos
# Row-standardization 
w.transform = "r"

# Corrigir para ilhas
#w.set_transform('r')
islands = w.islands
if islands:
    for island in islands:
        w.neighbors[island] = [island]
        w.weights[island] = [1]  
        
# Selecionar a coluna com a percentagem de população superior a 65 anos
data = gdf1315['IND65']


# Calcular a estatística de Moran
moran = esda.Moran(data, w)

# Imprimir o valor I de Moran e o valor-p
print('I de Moran:', moran.I)
print('Valor-p:', moran.p_sim)

      BGRI2021 DTMNFR21  N_INDIVIDUOS_65_OU_MAIS  N_INDIVIDUOS     IND65  \
0  13150600202   131506                      3.0          39.0  0.076923   
1  13150300506   131503                      7.0          36.0  0.194444   
2  13150600203   131506                      9.0          38.0  0.236842   
3  13150600204   131506                      3.0          14.0  0.214286   
4  13150300907   131503                     87.0         497.0  0.175050   
5  13150301501   131503                     11.0          24.0  0.458333   
6  13150302404   131503                     60.0         262.0  0.229008   
7  13150300601   131503                     11.0          51.0  0.215686   
8  13150301601   131503                     37.0         125.0  0.296000   
9  13150301102   131503                     44.0         188.0  0.234043   

      IND14     IND_H     IND_M  EDIF_3PISOS  
0  0.076923  0.487179  0.512821     0.083333  
1  0.111111  0.500000  0.500000     0.000000  
2  0.131579  0.526316  0.473684     0.066667  
3  0.142857  0.500000  0.500000     0.666667  
4  0.152918  0.496982  0.503018     0.067485  
5  0.083333  0.458333  0.541667     0.125000  
6  0.125954  0.442748  0.557252     0.454545  
7  0.098039  0.529412  0.470588     0.444444  
8  0.056000  0.456000  0.544000     0.500000  
9  0.117021  0.494681  0.505319     0.812500  

I de Moran: nan
Valor-p: 0.001

# Calcular o Spatial Lag
# 
gdf1315["w_IND14"] = weights.lag_spatial(w, gdf1315['IND14'])

# And their respective centered versions, where we subtract the average off of every value
gdf1315["IND14_std"] = gdf1315["IND14"] - gdf1315["IND14"].mean()
gdf1315["w_IND14_std"] = weights.lag_spatial(w, gdf1315['IND14_std'])

# Visualizar Valores em Cima e Baixo do Médio
# Set up the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values

sns.regplot(
    x="IND14_std", y="w_IND14_std", data=gdf1315 ci=None
)


# Criar os Quadrantos
# Set up the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values
sns.regplot(
    x="IND14_std", y="w_IND14_std", data=gdf1315, ci=None
)

# Esta Parte demora muito tempo

# Add vertical and horizontal lines (definição valor onde adicionar)
plt.axvline(0, c="k", alpha=0.5)
plt.axhline(0, c="k", alpha=0.5)

# # Adicionar Text para Cada Quadrant - coordinados Text tendo em conta a distribuição dos dados
# plt.text(0.6, 0.20, "HH", fontsize=25, c="r")
# plt.text(0.6, -0.15, "HL", fontsize=25, c="r")
# plt.text(-0.2, 0.20, "LH", fontsize=25, c="r")
# plt.text(-0.15, -0.15, "LL", fontsize=25, c="r")

# Displaby
plt.show()

9.4 Clustering e Regionalization

import matplotlib.pyplot as plt  # Graphics
from matplotlib import colors
import seaborn as sns  # Graphics
import geopandas as gpd # Spatial data manipulation
import pandas as pd  # Tabular data manipulation
import numpy as np

# Para Evitar Aviso Point Patterns
from shapely.geometry import Point
import contextily  # Background tiles

# Bibliotecas Referidos no Notebook
from esda.moran import Moran
from libpysal.weights import Queen, KNN
#import pysal.lib # importação geral
from pysal.explore import esda  # Exploratory Spatial analytics
from pysal.lib import weights  # Spatial weights

preparação dos dados

import geopandas as gpd
import matplotlib.pyplot as plt

# Caminho para o arquivo GeoPackage
gpk = r'data\geo\BGRI2021_1106.gpkg'

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1106 = gpd.read_file(gpk)

# Simplificar a geografia para uma precisão de 5 metros
gdf1106['geometry'] = gdf1106['geometry'].simplify(tolerance=5)

# Visualizar o GeoDataFrame
gdf1106.plot(column = 'DTMNFR21',
              legend = False)
              
# Nomes Atributos
#gdf1106.info()

<Axes: >

# Incluir os nomes das variáveis que devem ser utilizados para criar os agrupamentos (clusters)
cluster_variables = [
    gdf1106.columns[13],  # 
    gdf1106.columns[14],  # 
    gdf1106.columns[15],  # 
    gdf1106.columns[31],  # 
    gdf1106.columns[32],  # 
    gdf1106.columns[37],  # 
    gdf1106.columns[41],  # 
    gdf1106.columns[42],  # E
    gdf1106.columns[44] # 
]

# Tratamento NaN
# Contar o número total de NaNs no DataFrame
total_nans = gdf1106.isna().sum().sum()
print('Número total de registros com NaN:', total_nans)

# Contar o número de NaNs em cada coluna
nans_por_coluna = gdf1106.isna().sum()
print('Número de registros com NaN por coluna:\n', nans_por_coluna)

# Preencher com valor 0
gdf1106 = gdf1106.fillna(0)

Número total de registros com NaN: 0
Número de registros com NaN por coluna:
 OBJECTID                                                         0
BGRI2021                                                         0
DT21                                                             0
DTMN21                                                           0
DTMNFR21                                                         0
DTMNFRSEC21                                                      0
SECNUM21                                                         0
SSNUM21                                                          0
SECSSNUM21                                                       0
SUBSECCAO                                                        0
NUTS1                                                            0
NUTS2                                                            0
NUTS3                                                            0
N_EDIFICIOS_CLASSICOS                                            0
N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ                              0
N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS                    0
N_EDIFICIOS_EXCLUSIV_RESID                                       0
N_EDIFICIOS_1_OU_2_PISOS                                         0
N_EDIFICIOS_3_OU_MAIS_PISOS                                      0
N_EDIFICIOS_CONSTR_ANTES_1945                                    0
N_EDIFICIOS_CONSTR_1946_1980                                     0
N_EDIFICIOS_CONSTR_1981_2000                                     0
N_EDIFICIOS_CONSTR_2001_2010                                     0
N_EDIFICIOS_CONSTR_2011_2021                                     0
N_EDIFICIOS_COM_NECESSIDADES_REPARACAO                           0
N_ALOJAMENTOS_TOTAL                                              0
N_ALOJAMENTOS_FAMILIARES                                         0
N_ALOJAMENTOS_FAM_CLASS_RHABITUAL                                0
N_ALOJAMENTOS_FAM_CLASS_VAGOS_OU_RESID_SECUNDARIA                0
N_RHABITUAL_ACESSIVEL_CADEIRAS_RODAS                             0
N_RHABITUAL_COM_ESTACIONAMENTO                                   0
N_RHABITUAL_PROP_OCUP                                            0
N_RHABITUAL_ARRENDADOS                                           0
N_AGREGADOS_DOMESTICOS_PRIVADOS                                  0
N_ADP_1_OU_2_PESSOAS                                             0
N_ADP_3_OU_MAIS_PESSOAS                                          0
N_NUCLEOS_FAMILIARES                                             0
N_NUCLEOS_FAMILIARES_COM_FILHOS_TENDO_O_MAIS_NOVO_MENOS_DE_25    0
N_INDIVIDUOS                                                     0
N_INDIVIDUOS_H                                                   0
N_INDIVIDUOS_M                                                   0
N_INDIVIDUOS_0_14                                                0
N_INDIVIDUOS_15_24                                               0
N_INDIVIDUOS_25_64                                               0
N_INDIVIDUOS_65_OU_MAIS                                          0
SHAPE_Length                                                     0
SHAPE_Area                                                       0
geometry                                                         0
dtype: int64

# Mostrar como mapas temáticos os valores dos atributos escolhidos  

f, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 12))
# Make the axes accessible with single indexing
axs = axs.flatten()
# Start a loop over all the variables of interest
for i, col in enumerate(cluster_variables):
    # select the axis where the map will go
    ax = axs[i]
    # Plot the map
    gdf1106.plot(
        column=col,
        ax=ax,
        scheme="Quantiles",
        linewidth=0,
        cmap="RdPu",
    )
    # Remove axis clutter
    ax.set_axis_off()
    # Set the axis title to the name of variable being plotted
    ax.set_title(col)
# Display the figure
plt.show()

C:\Users\bruno.lima\AppData\Roaming\Python\Python312\site-packages\mapclassify\classifiers.py:1592: UserWarning: Not enough unique values in array to form 5 classes. Setting k to 4.
  self.bins = quantile(y, k=k)

Calcular Moran I para todas as variáveis

# Generate W from the GeoDataFrame
# w = weights.distance.KNN.from_dataframe(gdf1106, k=8)
# Metodo Alternativo:
w = weights.Queen.from_dataframe(gdf1106, use_index=True)

# Set seed for reproducibility
np.random.seed(123456)
# Calculate Moran's I for each variable
mi_results = [
    Moran(gdf1106[variable], w) for variable in cluster_variables
    ]
# Structure results as a list of tuples
mi_results = [
    (variable, res.I, res.p_sim)
    for variable, res in zip(cluster_variables, mi_results)
]
# Display on table
table = pd.DataFrame(
    mi_results, columns=["Variable", "Moran's I", "P-value"]
).set_index("Variable")
table

	Moran's I	P-value
Variable
N_EDIFICIOS_CLASSICOS	0.300793	0.001
N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ	0.384476	0.001
N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS	0.434427	0.001
N_RHABITUAL_PROP_OCUP	0.420872	0.001
N_RHABITUAL_ARRENDADOS	0.370622	0.001
N_NUCLEOS_FAMILIARES_COM_FILHOS_TENDO_O_MAIS_NOVO_MENOS_DE_25	0.391527	0.001
N_INDIVIDUOS_0_14	0.359343	0.001
N_INDIVIDUOS_15_24	0.371393	0.001
N_INDIVIDUOS_65_OU_MAIS	0.420039	0.001

# Utilizar kdeplot dá um aviso
import warnings

# Nao mostrar aviso FutereWarning (não aconselhável)
warnings.filterwarnings("ignore", category=FutureWarning)

# Mostrar kdeplot para cada variável
# _: Convenção Python que mostra que o resultado não está a ser utilizado

sns.pairplot(
    gdf1106[cluster_variables], kind="reg", diag_kind="kde"
)

# The distance between observations in terms of these variates can be computed easily using
from sklearn import metrics
metrics.pairwise_distances(
    gdf1106[[cluster_variables[0], cluster_variables[5]]].head()
).round(4)

array([[ 0.    ,  8.6023, 13.0384, 13.8924, 12.2066],
       [ 8.6023,  0.    , 18.1108, 19.105 , 17.1172],
       [13.0384, 18.1108,  0.    ,  1.    ,  1.    ],
       [13.8924, 19.105 ,  1.    ,  0.    ,  2.    ],
       [12.2066, 17.1172,  1.    ,  2.    ,  0.    ]])

from sklearn.preprocessing import robust_scale
# And create the db_scaled object which contains only the variables we are interested in, scaled:
db_scaled = robust_scale(gdf1106[cluster_variables])
print(type(db_scaled))

<class 'numpy.ndarray'>

9.4.1 Cluster GeoDemograficos

9.4.2 K-means

# Initialize KMeans instance
from sklearn.cluster import KMeans

# Initialize KMeans instance
# Definir n_init explicitemente
kmeans = KMeans(n_clusters=5,n_init=10)

# Set the seed for reproducibility
np.random.seed(1234)
# Run K-Means algorithm
k5cls = kmeans.fit(db_scaled)

# Print first five labels
k5cls.labels_[:5]

array([3, 3, 1, 1, 1])

visualização dos resultados

# Assign labels into a column
gdf1106["k5cls"] = k5cls.labels_
# Set up figure and ax
f, ax = plt.subplots(1, figsize=(9, 9))
# Plot unique values choropleth including
# a legend and with no boundary lines
gdf1106.plot(
    column="k5cls", categorical=True, legend=True, linewidth=0, ax=ax,
    legend_kwds={'loc': 'upper left'}
)
# Remove axis
ax.set_axis_off()
# Display the map
plt.show()

Análise dos resultados obtidos

# Group data table by cluster label and count observations
k5sizes = gdf1106.groupby("k5cls").size()
k5sizes

k5cls
0     269
1    1225
2     110
3     841
4     377
dtype: int64

obter a área de cada grupo (dissolve) utilizando o atributo SHAPE_Area

# Dissolve areas by Cluster, aggregate by summing,
# and keep column for area
areas = gdf1106.dissolve(by="k5cls", aggfunc="sum")["SHAPE_Area"]
areas

k5cls
0    1.147251e+07
1    5.281033e+07
2    7.158524e+06
3    1.926728e+07
4    9.345682e+06
Name: SHAPE_Area, dtype: float64

nº de elementos vs área de cada grupo

# Visualizar o nº de elementos vs a Área total

# Bind cluster figures in a single table
area_tracts = pd.DataFrame({"No. Tracts": k5sizes, "Area": areas})
# Convert raw values into percentages
area_tracts = area_tracts * 100 / area_tracts.sum()
# Bar plot
ax = area_tracts.plot.bar()
# Rename axes
ax.set_xlabel("Cluster labels")
ax.set_ylabel("Percentage by cluster")

plt.show()

Criar um perfil de cada cluster, executamos as seguintes tarefas:

1. Calcular os médios de cada variável
2. Arrumar (tidy up) os dados, assegurando que cada linha é uma observação e cada coluna é uma variável
   
3. Visualizar a distribuição dos valores das variáveis por grupo

# Criar os Perfis de cada Cluster em relação as variáveis de input
# Mostrar os médios das variáveis em cada cluster
# Group table by cluster label, keep the variables used
# for clustering, and obtain their mean
k5means = gdf1106.groupby("k5cls")[cluster_variables].mean()
# Transpose the table and print it rounding each value
# to three decimals
k5means.T.round(3)

k5cls	0	1	2	3	4
N_EDIFICIOS_CLASSICOS	27.112	7.980	63.055	17.776	27.236
N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ	2.584	2.332	51.491	1.860	20.599
N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS	24.227	5.529	11.318	15.744	6.552
N_RHABITUAL_PROP_OCUP	138.301	13.353	41.400	66.421	20.939
N_RHABITUAL_ARRENDADOS	114.539	13.180	49.818	50.599	19.419
N_NUCLEOS_FAMILIARES_COM_FILHOS_TENDO_O_MAIS_NOVO_MENOS_DE_25	85.810	7.092	28.118	36.056	12.011
N_INDIVIDUOS_0_14	86.721	7.330	28.955	36.922	12.406
N_INDIVIDUOS_15_24	67.617	5.908	21.355	28.056	9.812
N_INDIVIDUOS_65_OU_MAIS	140.792	15.361	53.645	66.652	24.199

# Index db on cluster ID
tidy_db = gdf1106.set_index("k5cls")
# Keep only variables used for clustering
tidy_db = tidy_db[cluster_variables]
# Stack column names into a column, obtaining
# a "long" version of the dataset
tidy_db = tidy_db.stack()
# Take indices into proper columns
tidy_db = tidy_db.reset_index()
# Rename column names
tidy_db = tidy_db.rename(
    columns={"level_1": "Attribute", 0: "Values"}
)
# Check out result
tidy_db.head()

	k5cls	Attribute	Values
0	3	N_EDIFICIOS_CLASSICOS	15.0
1	3	N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ	1.0
2	3	N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS	14.0
3	3	N_RHABITUAL_PROP_OCUP	30.0
4	3	N_RHABITUAL_ARRENDADOS	45.0

# hows the distribution of each cluster’s values for each variable. 
# This gives us the full distributional profile of each cluster:
# Scale fonts to make them more readable
sns.set(font_scale=1.5)
# Setup the facets
facets = sns.FacetGrid(
    data=tidy_db,
    col="Attribute",
    hue="k5cls",
    sharey=False,
    sharex=False,
    aspect=2,
    col_wrap=3,
)
# Build the plot from `sns.kdeplot`
facets.map(sns.kdeplot, "Values", shade=True).add_legend()

9.4.3 Hierarchical Clustering

from sklearn.cluster import AgglomerativeClustering

# Set seed for reproducibility
np.random.seed(0)
# Initialize the algorithm
model = AgglomerativeClustering(linkage="ward", n_clusters=5)
# Run clustering (input dataset é sempre o db_scaled com valores estandarizados)
model.fit(db_scaled)
# Assign labels to main data table
gdf1106["ward5"] = model.labels_

ward5sizes = gdf1106.groupby("ward5").size()
ward5sizes

ward5
0     578
1      93
2    1340
3     151
4     660
dtype: int64

visualizar dados

# Index db on cluster ID
tidy_db = gdf1106.set_index("ward5")
# Keep only variables used for clustering
tidy_db = tidy_db[cluster_variables]
# Stack column names into a column, obtaining
# a "long" version of the dataset
tidy_db = tidy_db.stack()
# Take indices into proper columns
tidy_db = tidy_db.reset_index()
# Rename column names
tidy_db = tidy_db.rename(
    columns={"level_1": "Attribute", 0: "Values"}
)
# Check out result
tidy_db.head()

	ward5	Attribute	Values
0	2	N_EDIFICIOS_CLASSICOS	15.0
1	2	N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ	1.0
2	2	N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS	14.0
3	2	N_RHABITUAL_PROP_OCUP	30.0
4	2	N_RHABITUAL_ARRENDADOS	45.0

# Setup the facets
facets = sns.FacetGrid(
    data=tidy_db,
    col="Attribute",
    hue="ward5",
    sharey=False,
    sharex=False,
    aspect=2,
    col_wrap=3,
)
# Build the plot as a `sns.kdeplot`
facets.map(sns.kdeplot, "Values", shade=True).add_legend()

comparação dos 2 resultados

gdf1106["ward5"] = model.labels_
# Set up figure and ax
f, axs = plt.subplots(1, 2, figsize=(12, 6))

### K-Means ###
ax = axs[0]
# Plot unique values choropleth including
# a legend and with no boundary lines
gdf1106.plot(
    column="ward5",
    categorical=True,
    cmap="Set3",
    legend=True,
    linewidth=0,
    ax=ax,
    legend_kwds={'loc': 'upper left'}
)
# Remove axis
ax.set_axis_off()
# Add title
ax.set_title("K-Means solution ($k=5$)")

### AHC ###
ax = axs[1]
# Plot unique values choropleth including
# a legend and with no boundary lines
gdf1106.plot(
    column="k5cls",
    categorical=True,
    cmap="Set3",
    legend=True,
    linewidth=0,
    ax=ax,
    legend_kwds={'loc': 'upper left'}
)
# Remove axis
ax.set_axis_off()
# Add title
ax.set_title("AHC solution ($k=5$)")

# Display the map
plt.show()

9.4.4 Regionalization

Criar novos Cluster

• Incluir matriz de vizinhança
• connectivity = w.sparse

# Set the seed for reproducibility
np.random.seed(123456)
# Specify cluster model with spatial constraint
model = AgglomerativeClustering(
    linkage="ward", connectivity=w.sparse, n_clusters=5
)
# Fit algorithm to the data
model.fit(db_scaled)

AgglomerativeClustering(connectivity=<2822x2822 sparse matrix of type '<class 'numpy.float64'>'
    with 15732 stored elements in Compressed Sparse Row format>,
                        n_clusters=5)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

visualizar resultados

gdf1106["ward5wq"] = model.labels_
# Set up figure and ax
f, ax = plt.subplots(1, figsize=(9, 9))
# Plot unique values choropleth including a legend and with no boundary lines
gdf1106.plot(
    column="ward5wq",
    categorical=True,
    legend=True,
    linewidth=0,
    ax=ax,
    legend_kwds={'loc': 'upper left'}
)
# Remove axis
ax.set_axis_off()
# Display the map
plt.show()

# Guardar o mapa:
f.savefig(r"data\geo\mapa_clusters1.png")

Repetir o processo com outra matriz de vizinhança

# Generate W from the GeoDataFrame
w = weights.distance.KNN.from_dataframe(gdf1106, k=8)
# Metodo Alternativo:
#w = weights.Queen.from_dataframe(gdf1106, use_index=True)

9.4.5 Outros metodos

Medida compactness

results = []
for cluster_type in ("k5cls", "ward5", "ward5wq", "ward5wknn"):
    # compute the region polygons using a dissolve
    regions = db[[cluster_type, "geometry"]].dissolve(by=cluster_type)
    # compute the actual isoperimetric quotient for these regions
    ipqs = (
        regions.area * 4 * numpy.pi / (regions.boundary.length ** 2)
    )
    # cast to a dataframe
    result = ipqs.to_frame(cluster_type)
    results.append(result)
# stack the series together along columns
pandas.concat(results, axis=1)

import numpy
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

# gdf1106['IND65'] = gdf1106.N_INDIVIDUOS_65_OU_MAIS/gdf1106.N_INDIVIDUOS
gdf1106['compact'] = gdf1106.geometry.area * 4  * numpy.pi / (gdf1106.geometry.boundary.length ** 2)

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 3}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf1106.plot(column='compact', 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=9, 
                      cmap='PuBu', 
                      legend=True,
                      edgecolor = 'None', # sem outline
                      legend_kwds  = lgnd_kwds)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Compactness')
plt.show()

print(gdf1106.head())

   OBJECTID     BGRI2021 DT21 DTMN21 DTMNFR21 DTMNFRSEC21 SECNUM21 SSNUM21  \
0     30243  11065602301   11   1106   110656   110656023      023      01   
1     30244  11065700203   11   1106   110657   110657002      002      03   
2     30311  11065801011   11   1106   110658   110658010      010      11   
3     30312  11065801012   11   1106   110658   110658010      010      12   
4     30313  11065801013   11   1106   110658   110658010      010      13   

  SECSSNUM21    SUBSECCAO  ... N_INDIVIDUOS_15_24 N_INDIVIDUOS_25_64  \
0      02301  11065602301  ...               31.0              127.0   
1      00203  11065700203  ...               24.0               85.0   
2      01011  11065801011  ...                5.0                9.0   
3      01012  11065801012  ...                2.0                7.0   
4      01013  11065801013  ...                3.0                8.0   

  N_INDIVIDUOS_65_OU_MAIS  SHAPE_Length   SHAPE_Area  \
0                    30.0    409.853268  9657.766943   
1                    62.0    368.048569  8445.226728   
2                     3.0    239.674582  2582.412127   
3                     7.0    250.811945  3077.844919   
4                     4.0    238.956181  2555.620562   

                                            geometry  k5cls  ward5  ward5wq  \
0  POLYGON ((-86809.545 -103264.238, -86801.039 -...      3      2        2   
1  POLYGON ((-88183.921 -103236.850, -88218.744 -...      3      2        2   
2  POLYGON ((-94424.359 -107038.246, -94495.548 -...      1      2        0   
3  POLYGON ((-94547.783 -107023.019, -94626.237 -...      1      2        0   
4  POLYGON ((-94440.784 -107017.804, -94514.143 -...      1      2        0   

    compact  
0  0.729871  
1  0.783449  
2  0.564926  
3  0.611801  
4  0.562432  

[5 rows x 52 columns]

Kriging com o package pykrige

import numpy as np
from pykrige.ok import OrdinaryKriging
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt  

# Sample data points
data = np.array(
    [
        [0.3, 1.2, 0.5],
        [1.1, 3.2, 0.4],
        [1.8, 0.8, 0.6],
        [2.8, 2.6, 0.7],
        [3.2, 0.3, 0.8],
    ]
)  # [x, y, z]

# Define the grid to interpolate onto
gridx = np.arange(0.0, 4.1, 0.1)
gridy = np.arange(0.0, 4.1, 0.1)

# Create an Ordinary Kriging object
OK = OrdinaryKriging(
    data[:, 0],  # X coordinates
    data[:, 1],  # Y coordinates
    data[:, 2],  # Z values
    variogram_model="spherical",  # Variogram model (can also use "linear" gaussian" or "spherical")
    verbose=False,
    enable_plotting=True,  # Enable plotting of the variogram (optional)
)

# Execute Ordinary Kriging on the defined grid
# `z` contains the interpolated values
# `ss` contains the standard deviation at each grid point
z, ss = OK.execute("grid", gridx, gridy)

# Writes the kriged grid to an ASCII grid file and plot it.
kt.write_asc_grid(gridx, gridy, z, filename="data\geo\output.asc")
plt.imshow(z)
plt.show()

<>:37: SyntaxWarning: invalid escape sequence '\g'
<>:37: SyntaxWarning: invalid escape sequence '\g'
C:\Users\bruno.lima\AppData\Local\Temp\ipykernel_9296\373054154.py:37: SyntaxWarning: invalid escape sequence '\g'
  kt.write_asc_grid(gridx, gridy, z, filename="data\geo\output.asc")

9.4.6 Exercício

adaptar o código anterior para outro municipio

preparação dos dados

import geopandas as gpd
import matplotlib.pyplot as plt

# Caminho para o arquivo GeoPackage municipio Valongo
gpk = r'data\geo\BGRI2021_1315.gpkg'

# Ler os dados do GeoPackage para um GeoDataFrame
gdf1315 = gpd.read_file(gpk)

# Simplificar a geografia para uma precisão de 5 metros
gdf1315['geometry'] = gdf1315['geometry'].simplify(tolerance=5)

# # Visualizar o GeoDataFrame
# gdf1315.plot(column = 'DTMNFR21',
#               legend = False)
              
# Nomes Atributos
gdf1315.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 1001 entries, 0 to 1000
Data columns (total 48 columns):
 #   Column                                                         Non-Null Count  Dtype   
---  ------                                                         --------------  -----   
 0   OBJECTID                                                       1001 non-null   int64   
 1   BGRI2021                                                       1001 non-null   object  
 2   DT21                                                           1001 non-null   object  
 3   DTMN21                                                         1001 non-null   object  
 4   DTMNFR21                                                       1001 non-null   object  
 5   DTMNFRSEC21                                                    1001 non-null   object  
 6   SECNUM21                                                       1001 non-null   object  
 7   SSNUM21                                                        1001 non-null   object  
 8   SECSSNUM21                                                     1001 non-null   object  
 9   SUBSECCAO                                                      1001 non-null   object  
 10  NUTS1                                                          1001 non-null   object  
 11  NUTS2                                                          1001 non-null   object  
 12  NUTS3                                                          1001 non-null   object  
 13  N_EDIFICIOS_CLASSICOS                                          1001 non-null   float64 
 14  N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ                            1001 non-null   float64 
 15  N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS                  1001 non-null   float64 
 16  N_EDIFICIOS_EXCLUSIV_RESID                                     1001 non-null   float64 
 17  N_EDIFICIOS_1_OU_2_PISOS                                       1001 non-null   float64 
 18  N_EDIFICIOS_3_OU_MAIS_PISOS                                    1001 non-null   float64 
 19  N_EDIFICIOS_CONSTR_ANTES_1945                                  1001 non-null   float64 
 20  N_EDIFICIOS_CONSTR_1946_1980                                   1001 non-null   float64 
 21  N_EDIFICIOS_CONSTR_1981_2000                                   1001 non-null   float64 
 22  N_EDIFICIOS_CONSTR_2001_2010                                   1001 non-null   float64 
 23  N_EDIFICIOS_CONSTR_2011_2021                                   1001 non-null   float64 
 24  N_EDIFICIOS_COM_NECESSIDADES_REPARACAO                         1001 non-null   float64 
 25  N_ALOJAMENTOS_TOTAL                                            1001 non-null   float64 
 26  N_ALOJAMENTOS_FAMILIARES                                       1001 non-null   float64 
 27  N_ALOJAMENTOS_FAM_CLASS_RHABITUAL                              1001 non-null   float64 
 28  N_ALOJAMENTOS_FAM_CLASS_VAGOS_OU_RESID_SECUNDARIA              1001 non-null   float64 
 29  N_RHABITUAL_ACESSIVEL_CADEIRAS_RODAS                           1001 non-null   float64 
 30  N_RHABITUAL_COM_ESTACIONAMENTO                                 1001 non-null   float64 
 31  N_RHABITUAL_PROP_OCUP                                          1001 non-null   float64 
 32  N_RHABITUAL_ARRENDADOS                                         1001 non-null   float64 
 33  N_AGREGADOS_DOMESTICOS_PRIVADOS                                1001 non-null   float64 
 34  N_ADP_1_OU_2_PESSOAS                                           1001 non-null   float64 
 35  N_ADP_3_OU_MAIS_PESSOAS                                        1001 non-null   float64 
 36  N_NUCLEOS_FAMILIARES                                           1001 non-null   float64 
 37  N_NUCLEOS_FAMILIARES_COM_FILHOS_TENDO_O_MAIS_NOVO_MENOS_DE_25  1001 non-null   float64 
 38  N_INDIVIDUOS                                                   1001 non-null   float64 
 39  N_INDIVIDUOS_H                                                 1001 non-null   float64 
 40  N_INDIVIDUOS_M                                                 1001 non-null   float64 
 41  N_INDIVIDUOS_0_14                                              1001 non-null   float64 
 42  N_INDIVIDUOS_15_24                                             1001 non-null   float64 
 43  N_INDIVIDUOS_25_64                                             1001 non-null   float64 
 44  N_INDIVIDUOS_65_OU_MAIS                                        1001 non-null   float64 
 45  SHAPE_Length                                                   1001 non-null   float64 
 46  SHAPE_Area                                                     1001 non-null   float64 
 47  geometry                                                       1001 non-null   geometry
dtypes: float64(34), geometry(1), int64(1), object(12)
memory usage: 375.5+ KB

# Incluir os nomes das variáveis que devem ser utilizados para criar os agrupamentos (clusters)
cluster_variables = [
    gdf1315.columns[24],  # 
    gdf1315.columns[25],  # 
    gdf1315.columns[26],  # 
    gdf1315.columns[27],  # 
    gdf1315.columns[28],
    gdf1315.columns[38],
    gdf1315.columns[39],
    gdf1315.columns[40]
]

# Tratamento NaN
# Contar o número total de NaNs no DataFrame
total_nans = gdf1315.isna().sum().sum()
print('Número total de registros com NaN:', total_nans)

# Contar o número de NaNs em cada coluna
nans_por_coluna = gdf1315.isna().sum()
print('Número de registros com NaN por coluna:\n', nans_por_coluna)

# Preencher com valor 0
gdf1315 = gdf1315.fillna(0)

Número total de registros com NaN: 0
Número de registros com NaN por coluna:
 OBJECTID                                                         0
BGRI2021                                                         0
DT21                                                             0
DTMN21                                                           0
DTMNFR21                                                         0
DTMNFRSEC21                                                      0
SECNUM21                                                         0
SSNUM21                                                          0
SECSSNUM21                                                       0
SUBSECCAO                                                        0
NUTS1                                                            0
NUTS2                                                            0
NUTS3                                                            0
N_EDIFICIOS_CLASSICOS                                            0
N_EDIFICIOS_CLASS_CONST_1_OU_2_ALOJ                              0
N_EDIFICIOS_CLASS_CONST_3_OU_MAIS_ALOJAMENTOS                    0
N_EDIFICIOS_EXCLUSIV_RESID                                       0
N_EDIFICIOS_1_OU_2_PISOS                                         0
N_EDIFICIOS_3_OU_MAIS_PISOS                                      0
N_EDIFICIOS_CONSTR_ANTES_1945                                    0
N_EDIFICIOS_CONSTR_1946_1980                                     0
N_EDIFICIOS_CONSTR_1981_2000                                     0
N_EDIFICIOS_CONSTR_2001_2010                                     0
N_EDIFICIOS_CONSTR_2011_2021                                     0
N_EDIFICIOS_COM_NECESSIDADES_REPARACAO                           0
N_ALOJAMENTOS_TOTAL                                              0
N_ALOJAMENTOS_FAMILIARES                                         0
N_ALOJAMENTOS_FAM_CLASS_RHABITUAL                                0
N_ALOJAMENTOS_FAM_CLASS_VAGOS_OU_RESID_SECUNDARIA                0
N_RHABITUAL_ACESSIVEL_CADEIRAS_RODAS                             0
N_RHABITUAL_COM_ESTACIONAMENTO                                   0
N_RHABITUAL_PROP_OCUP                                            0
N_RHABITUAL_ARRENDADOS                                           0
N_AGREGADOS_DOMESTICOS_PRIVADOS                                  0
N_ADP_1_OU_2_PESSOAS                                             0
N_ADP_3_OU_MAIS_PESSOAS                                          0
N_NUCLEOS_FAMILIARES                                             0
N_NUCLEOS_FAMILIARES_COM_FILHOS_TENDO_O_MAIS_NOVO_MENOS_DE_25    0
N_INDIVIDUOS                                                     0
N_INDIVIDUOS_H                                                   0
N_INDIVIDUOS_M                                                   0
N_INDIVIDUOS_0_14                                                0
N_INDIVIDUOS_15_24                                               0
N_INDIVIDUOS_25_64                                               0
N_INDIVIDUOS_65_OU_MAIS                                          0
SHAPE_Length                                                     0
SHAPE_Area                                                       0
geometry                                                         0
dtype: int64

# Mostrar como mapas temáticos os valores dos atributos escolhidos  

f, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 12))
# Make the axes accessible with single indexing
axs = axs.flatten()
# Start a loop over all the variables of interest
for i, col in enumerate(cluster_variables):
    # select the axis where the map will go
    ax = axs[i]
    # Plot the map
    gdf1315.plot(
        column=col,
        ax=ax,
        scheme="Quantiles",
        linewidth=0,
        cmap="RdPu",
    )
    # Remove axis clutter
    ax.set_axis_off()
    # Set the axis title to the name of variable being plotted
    ax.set_title(col)
# Display the figure
plt.show()

Calcular Moran I para todas as variáveis

# Generate W from the GeoDataFrame
# w = weights.distance.KNN.from_dataframe(gdf1315, k=8)
# Metodo Alternativo:
w = weights.Queen.from_dataframe(gdf1315, use_index=True)

# Set seed for reproducibility
np.random.seed(123456)
# Calculate Moran's I for each variable
mi_results = [
    Moran(gdf1315[variable], w) for variable in cluster_variables
    ]
# Structure results as a list of tuples
mi_results = [
    (variable, res.I, res.p_sim)
    for variable, res in zip(cluster_variables, mi_results)
]
# Display on table
table = pd.DataFrame(
    mi_results, columns=["Variable", "Moran's I", "P-value"]
).set_index("Variable")
table

	Moran's I	P-value
Variable
N_EDIFICIOS_COM_NECESSIDADES_REPARACAO	0.197635	0.001
N_ALOJAMENTOS_TOTAL	0.278853	0.001
N_ALOJAMENTOS_FAMILIARES	0.278712	0.001
N_ALOJAMENTOS_FAM_CLASS_RHABITUAL	0.284710	0.001
N_ALOJAMENTOS_FAM_CLASS_VAGOS_OU_RESID_SECUNDARIA	0.167532	0.001
N_INDIVIDUOS	0.264845	0.001
N_INDIVIDUOS_H	0.255916	0.001
N_INDIVIDUOS_M	0.270961	0.001

# Utilizar kdeplot dá um aviso
import warnings

# Nao mostrar aviso FutereWarning (não aconselhável)
warnings.filterwarnings("ignore", category=FutureWarning)

# Mostrar kdeplot para cada variável
# _: Convenção Python que mostra que o resultado não está a ser utilizado

sns.pairplot(
    gdf1315[cluster_variables], kind="reg", diag_kind="kde"
)

# The distance between observations in terms of these variates can be computed easily using
from sklearn import metrics
metrics.pairwise_distances(
    gdf1315[[cluster_variables[0], cluster_variables[5]]].head()
).round(4)

array([[  0.    ,  12.3693,   9.0554,  25.    , 458.1572],
       [ 12.3693,   0.    ,   3.6056,  25.0599, 461.    ],
       [  9.0554,   3.6056,   0.    ,  25.632 , 459.0098],
       [ 25.    ,  25.0599,  25.632 ,   0.    , 483.149 ],
       [458.1572, 461.    , 459.0098, 483.149 ,   0.    ]])

from sklearn.preprocessing import robust_scale
# And create the db_scaled object which contains only the variables we are interested in, scaled:
db_scaled = robust_scale(gdf1315[cluster_variables])
print(type(db_scaled))

<class 'numpy.ndarray'>

import numpy
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

# gdf1106['IND65'] = gdf1106.N_INDIVIDUOS_65_OU_MAIS/gdf1106.N_INDIVIDUOS
gdf1315['compact'] = gdf1315.geometry.area * 4  * numpy.pi / (gdf1315.geometry.boundary.length ** 2)

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 3}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf1315.plot(column='compact', 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=9, 
                      cmap='PuBu', 
                      legend=True,
                      edgecolor = 'None', # sem outline
                      legend_kwds  = lgnd_kwds)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Compactness')
plt.show()

print(gdf1315.head())

   OBJECTID     BGRI2021 DT21 DTMN21 DTMNFR21 DTMNFRSEC21 SECNUM21 SSNUM21  \
0    125534  13150600202   13   1315   131506   131506002      002      02   
1    125535  13150300506   13   1315   131503   131503005      005      06   
2    125536  13150600203   13   1315   131506   131506002      002      03   
3    125537  13150600204   13   1315   131506   131506002      002      04   
4    125538  13150300907   13   1315   131503   131503009      009      07   

  SECSSNUM21    SUBSECCAO  ... N_INDIVIDUOS_H N_INDIVIDUOS_M  \
0      00202  13150600202  ...           19.0           20.0   
1      00506  13150300506  ...           18.0           18.0   
2      00203  13150600203  ...           20.0           18.0   
3      00204  13150600204  ...            7.0            7.0   
4      00907  13150300907  ...          247.0          250.0   

  N_INDIVIDUOS_0_14  N_INDIVIDUOS_15_24  N_INDIVIDUOS_25_64  \
0               3.0                 5.0                28.0   
1               4.0                 4.0                21.0   
2               5.0                 1.0                23.0   
3               2.0                 3.0                 6.0   
4              76.0                55.0               279.0   

   N_INDIVIDUOS_65_OU_MAIS  SHAPE_Length     SHAPE_Area  \
0                      3.0    523.018918    9289.852833   
1                      7.0   1109.669756   50141.267236   
2                      9.0    511.977438   15032.503593   
3                      3.0    553.145862   10797.658643   
4                     87.0   2328.672435  260550.923540   

                                            geometry   compact  
0  POLYGON ((-27481.127 171972.950, -27461.897 17...  0.411844  
1  POLYGON ((-34607.059 171934.156, -34553.823 17...  0.517489  
2  POLYGON ((-28357.470 171987.554, -28332.531 17...  0.715822  
3  POLYGON ((-27497.017 172000.233, -27486.203 17...  0.451444  
4  POLYGON ((-33995.356 171989.441, -33899.639 17...  0.611464  

[5 rows x 49 columns]

Kriging com o package pykrige

import numpy as np
from pykrige.ok import OrdinaryKriging
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt  

# Sample data points
data = np.array(
    [
        [0.3, 1.2, 0.5],
        [1.1, 3.2, 0.4],
        [1.8, 0.8, 0.6],
        [2.8, 2.6, 0.7],
        [3.2, 0.3, 0.8],
    ]
)  # [x, y, z]

# Define the grid to interpolate onto
gridx = np.arange(0.0, 4.1, 0.1)
gridy = np.arange(0.0, 4.1, 0.1)

# Create an Ordinary Kriging object
OK = OrdinaryKriging(
    data[:, 0],  # X coordinates
    data[:, 1],  # Y coordinates
    data[:, 2],  # Z values
    variogram_model="spherical",  # Variogram model (can also use "linear" gaussian" or "spherical")
    verbose=False,
    enable_plotting=True,  # Enable plotting of the variogram (optional)
)

# Execute Ordinary Kriging on the defined grid
# `z` contains the interpolated values
# `ss` contains the standard deviation at each grid point
z, ss = OK.execute("grid", gridx, gridy)

# Writes the kriged grid to an ASCII grid file and plot it.
kt.write_asc_grid(gridx, gridy, z, filename="data\geo\output.asc")
plt.imshow(z)
plt.show()

<>:37: SyntaxWarning: invalid escape sequence '\g'
<>:37: SyntaxWarning: invalid escape sequence '\g'
C:\Users\bruno.lima\AppData\Local\Temp\ipykernel_9296\373054154.py:37: SyntaxWarning: invalid escape sequence '\g'
  kt.write_asc_grid(gridx, gridy, z, filename="data\geo\output.asc")

9.5 JSON / API INE

9.5.1 JavaScript Object Notation

• formato de dados leve e eficiente
• dados estruturados
• amplamente utilizado

útil para:

• recolha de dados
• preparação de dados
• análise de dados
• visualização de dados

JSON vs XML

Um objecto JSON é um conjunto de pares de chave e valor encerrados por chaves.

packages python para ler JSON: json, pandas, requests, jsonpath-ng, jsonschema

# Ler JSON de um Dicionario
import pprint # para visualizar de forma mais amigavel os dados
import json

# Criar um Dictionary (exemplo dados municipios)
municipios = {
    "Almada": {
        "populacao": 177400        
    },
    "Cascais": {
        "populacao": 214124        
    },
    "Seixal": {
        "populacao": 160000
    },
    "Entroncamento": {
        "populacao": 20141
    },
    "Cadaval": {
        "populacao": 13372
    },
    "Sintra": {
        "populacao": 385606
    }
}

print(type(municipios))

# Converter o Dictionary para uma string JSON
mn_json = json.dumps(municipios)
print(type(mn_json))

# Fazer Load do JSON
data_json = json.loads(mn_json)
print(type(data_json))

# Mostrar População Cascais
print ("População Cascais:",data_json["Cascais"]["populacao"])

# Imprimir o Type do Object Devolvido
print('Tipo Objecto:', type(data_json))

# Imprimir o objecto com PrettyPrinter
pp = pprint.PrettyPrinter(indent=2)

# Imprimir o dicionário
pp.pprint(data_json)

<class 'dict'>
<class 'str'>
<class 'dict'>
População Cascais: 214124
Tipo Objecto: <class 'dict'>
{ 'Almada': {'populacao': 177400},
  'Cadaval': {'populacao': 13372},
  'Cascais': {'populacao': 214124},
  'Entroncamento': {'populacao': 20141},
  'Seixal': {'populacao': 160000},
  'Sintra': {'populacao': 385606}}

# Criar JSOn a partir de uma Listagem:
# import json
# import pprint

# Exemplo cidades de Portugal
cidades = ["Lisboa", "Porto", "Vila Nova de Gaia", "Amadora", "Braga", "Funchal", "Coimbra", "Almada", "Setúbal", "Agualva-Cacém"]

# Converter a lista para uma string JSON
cidades_json = json.dumps(cidades)
print('Tipo Objecto após json.dumps:', type(cidades_json))

# Converter a string JSON de volta para uma lista
data_json = json.loads(cidades_json)


# Imprimir o Type do Object Devolvido
print('Tipo Objecto após json.loads:', type(data_json))

# Imprimir os dados carregados
# Criar um objeto PrettyPrinter
pp = pprint.PrettyPrinter(indent=4)

# Imprimir o dicionário
pp.pprint(data_json)

Tipo Objecto após json.dumps: <class 'str'>
Tipo Objecto após json.loads: <class 'list'>
[   'Lisboa',
    'Porto',
    'Vila Nova de Gaia',
    'Amadora',
    'Braga',
    'Funchal',
    'Coimbra',
    'Almada',
    'Setúbal',
    'Agualva-Cacém']

ler dados de um ficheiro

# exemplo de um ficheiro com os municipios de Madrid
import json
import pandas as pd
# Ler Dados de um Ficheiro no computador (Municipios de Provincia de Madrid)
jsonfile = r"data\geo\municipio_comunidad_madrid.json"

# Abrir Ficheiro e fazer Load
with open(jsonfile, 'r') as f:
    json_data = json.load(f)

print(type(json_data))
    
# Verificar Tipo de Dados Devolvido e mostrar informação
if isinstance(json_data, list):
    print("JSON object is a list.")
    if json_obj:
        print("Numero Registos:", len(json_data))
        print("Registo Exemplo:", json_data[0])
elif isinstance(json_data, dict):
    print("JSON object is a dictionary.")
    print("Keys do Object:", list(json_data.keys()))
else:
    print("Unknown JSON object type.")

novoelem = json_data["data"]  
print (type(novoelem))

df = pd.DataFrame(novoelem)
print(df.info())

<class 'dict'>
JSON object is a dictionary.
Keys do Object: ['data']
<class 'list'>
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 179 entries, 0 to 178
Data columns (total 7 columns):
 #   Column                Non-Null Count  Dtype  
---  ------                --------------  -----  
 0   municipio_codigo      179 non-null    object 
 1   densidad_por_km2      179 non-null    float64
 2   municipio_codigo_ine  179 non-null    object 
 3   nuts4_nombre          179 non-null    object 
 4   municipio_nombre      179 non-null    object 
 5   nuts4_codigo          179 non-null    object 
 6   superficie_km2        179 non-null    float64
dtypes: float64(2), object(5)
memory usage: 9.9+ KB
None

ler dados de um URL

import json
import requests

proxies = {
  'http': 'http://proxy.ine.pt:8080',
  'https': 'http://proxy.ine.pt:8080',
}


url = "https://datos.comunidad.madrid/catalogo/dataset/032474a0-bf11-4465-bb92-392052962866/resource/301aed82-339b-4005-ab20-06db41ee7017/download/municipio_comunidad_madrid.json"

# Make an HTTP GET request to fetch the JSON data from the URL
# 
response = requests.get(url)#, proxies=proxies)

# Verificar Resposta
# Respostas Possiveis
if response.status_code == 200:
    #Obter JSON response
    json_data = response.json()
else:
    print("Failed to fetch data. Status code:", response.status_code)
    exit()


# Verificar Tipo de Dados Devolvido e mostrar informação
if isinstance(json_data, list):
    print("JSON object is a list.")
    if json_data:
        print("Numero Registos:", len(json_data))
        print("Registo Exemplo:", json_data[0])
elif isinstance(json_data, dict):
    print("JSON object is a dictionary.")
    print("Keys do Object:", list(json_data.keys()))
else:
    print("Unknown JSON object type.")

JSON object is a dictionary.
Keys do Object: ['data']

9.5.2 converter para pandas DataFrame

ler dados de uma listagem normal

  import pandas as pd
# Ler Informacao JSOn Municipios:
df_mn = pd.read_json(r'data\geo\municipios.json')

print(df_mn.info())
print('Primeiro Municipio',df_mn['municipio'][0])

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5 entries, 0 to 4
Data columns (total 2 columns):
 #   Column     Non-Null Count  Dtype 
---  ------     --------------  ----- 
 0   municipio  5 non-null      object
 1   populacao  5 non-null      int64 
dtypes: int64(1), object(1)
memory usage: 212.0+ bytes
None
Primeiro Municipio Almada

nested data

import pandas as pd
# Ler Informacao JSOn Municipios:
df_frmn = pd.read_json(r"data\geo\municipiosfreguesias.json")
print(df_frmn.info())
print('Primeiro Municipio',df_frmn['municipio'][0])

# Ver o tipo do atributo freguesias - Series
print(type(df_frmn['freguesias']))

# Selecionar o primeiro municipio
df_frmn = df_frmn.loc[df_frmn['municipio'] == 'Almada']

# Mostrar os valores das freguesias
print(df_frmn['freguesias'])

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5 entries, 0 to 4
Data columns (total 3 columns):
 #   Column      Non-Null Count  Dtype 
---  ------      --------------  ----- 
 0   municipio   5 non-null      object
 1   populacao   5 non-null      int64 
 2   freguesias  5 non-null      object
dtypes: int64(1), object(2)
memory usage: 252.0+ bytes
None
Primeiro Municipio Almada
<class 'pandas.core.series.Series'>
0    [{'nome': 'Almada', 'populacao': 53194}, {'nom...
Name: freguesias, dtype: object

utilizar json_normalize()

# Utilizar a funcao json_normalize para criar um DF apenas das freguesias

# Load the JSON file
with open(r"data\geo\municipiosfreguesias.json", "r", encoding="utf-8") as f:
    municipalities = json.load(f)

# Convert the JSON into a DataFrame
df = pd.json_normalize(municipalities, record_path=['freguesias'])

# Print the DataFrame
print(df)
print(df.info())

                   nome  populacao
0                Almada      53194
1              Cacilhas      20540
2              Caparica      29222
3     Costa da Caparica      33477
4    Laranjeiro e Feijó      19910
5   Marinha da Caparica      10090
6               Pichela      10012
7               Sobreda      10084
8           Alcabideche      26613
9   Carcavelos e Parede      74302
10              Cascais      36744
11              Estoril      40926
12              Guincho      24759
13                Amora      41588
14     Gandra e Samouco      13069
15              Palmela      23698
16               Seixal     120464
17        Entroncamento      39343
18              Cadaval      12889
19   Ferreira do Zêzere      10890
20             Lourinhã      11729
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 21 entries, 0 to 20
Data columns (total 2 columns):
 #   Column     Non-Null Count  Dtype 
---  ------     --------------  ----- 
 0   nome       21 non-null     object
 1   populacao  21 non-null     int64 
dtypes: int64(1), object(1)
memory usage: 468.0+ bytes
None

criar dataframe utilizando um loop

# Alternativa para criar Informação das Freguesias
# Ler os dados
with open(r"data\geo\municipiosfreguesias.json", "r", encoding="utf-8") as f:
    municipalities = json.load(f)
    
# Criar Lista das linhas do INPUT
rows = []

# Percorrer os municipios
for municipality in municipalities:
    # Iterate over the freguesias
    for freguesia in municipality['freguesias']:
        # Criar novo registo
        row = {
            'municipio': municipality['municipio'],
            'freguesia': freguesia['nome'],
            'populacao': freguesia['populacao']
        }

        # Adicionar Registo a Lista
        rows.append(row)

# Criar DataFame
df = pd.DataFrame(rows)

# Print the DataFrame
print(df)

        municipio            freguesia  populacao
0          Almada               Almada      53194
1          Almada             Cacilhas      20540
2          Almada             Caparica      29222
3          Almada    Costa da Caparica      33477
4          Almada   Laranjeiro e Feijó      19910
5          Almada  Marinha da Caparica      10090
6          Almada              Pichela      10012
7          Almada              Sobreda      10084
8         Cascais          Alcabideche      26613
9         Cascais  Carcavelos e Parede      74302
10        Cascais              Cascais      36744
11        Cascais              Estoril      40926
12        Cascais              Guincho      24759
13         Seixal                Amora      41588
14         Seixal     Gandra e Samouco      13069
15         Seixal              Palmela      23698
16         Seixal               Seixal     120464
17  Entroncamento        Entroncamento      39343
18        Cadaval              Cadaval      12889
19        Cadaval   Ferreira do Zêzere      10890
20        Cadaval             Lourinhã      11729

9.5.3 Exercício

A partir do objecto pt_info

1. Converter o dicionário para uma string JSON
   

2. Converter a string JSON de volta para um dicionário

3. Obter o tipo do objecto devolvido

4. Imprimir as chaves

5. Converter para DataFrame

6. Fazer algumas pesquisas

   • População de Vila Nova de Gaia
   • Area de Portugal

# Dictionary com informcao Portugal:
# Também é Nested
pt_info = {
  "country": "Portugal",
  "capital": "Lisbon",
  "population": 10347892,
  "area": 92212,
  "language": "Portuguese",
  "currency": "Euro",
  "cities": [
    {
      "name": "Lisbon",
      "population": 504762,
      "region": "Lisboa"
    },
    {
      "name": "Porto",
      "population": 219419,
      "region": "Norte"
    },
    {
      "name": "Vila Nova de Gaia",
      "population": 301877,
      "region": "Norte"
    },
    {
      "name": "Matosinhos",
      "population": 174339,
      "region": "Norte"
    },
    {
      "name": "Almada",
      "population": 174033,
      "region": "Lisboa"
    }
  ]
}

print(type(pt_info))

<class 'dict'>

# Converter o Dictionary para uma string JSON
mn_json = json.dumps(pt_info)
print(type(mn_json))

# Fazer Load do JSON
data_json = json.loads(mn_json)
print(type(data_json))

# Imprimir o Type do Object Devolvido
print('Tipo Objecto:', type(data_json))

# Imprimir o objecto com PrettyPrinter
pp = pprint.PrettyPrinter(indent=2)

# Imprimir o dicionário
pp.pprint(data_json)

# Obter o Nome da primeira cidade 
print(data_json['cities'][0]['name'])

# população do Porto
print('População do Porto: ',data_json['cities']['name' == 'Porto']['population'])

# 4. Imprimir as chaves 
# Verificar Tipo de Dados Devolvido e mostrar informação
if isinstance(data_json, list):
    print("JSON object is a list.")
    if json_obj:
        print("Numero Registos:", len(data_json))
        print("Registo Exemplo:", data_json[0])
elif isinstance(data_json, dict):
    print("JSON object is a dictionary.")
    print("Keys do Object:", list(data_json.keys()))
else:
    print("Unknown JSON object type.")

# 5. Converter para DataFrame
df = pd.json_normalize(data_json, record_path=['cities'])
print(df)


# 6. Fazer algumas pesquisas
#  - População de Vila Nova de Gaia
#  - Area de Portugal
populacao_vila_nova_de_gaia = df[df['name'] == 'Vila Nova de Gaia']['population'].values[0]
print("População de Vila Nova de Gaia:", populacao_vila_nova_de_gaia)

area_portugal = data_json['area']
print("Área de Portugal:", area_portugal)

<class 'str'>
<class 'dict'>
Tipo Objecto: <class 'dict'>
{ 'area': 92212,
  'capital': 'Lisbon',
  'cities': [ {'name': 'Lisbon', 'population': 504762, 'region': 'Lisboa'},
              {'name': 'Porto', 'population': 219419, 'region': 'Norte'},
              { 'name': 'Vila Nova de Gaia',
                'population': 301877,
                'region': 'Norte'},
              {'name': 'Matosinhos', 'population': 174339, 'region': 'Norte'},
              {'name': 'Almada', 'population': 174033, 'region': 'Lisboa'}],
  'country': 'Portugal',
  'currency': 'Euro',
  'language': 'Portuguese',
  'population': 10347892}
Lisbon
População do Porto:  504762
JSON object is a dictionary.
Keys do Object: ['country', 'capital', 'population', 'area', 'language', 'currency', 'cities']
                name  population  region
0             Lisbon      504762  Lisboa
1              Porto      219419   Norte
2  Vila Nova de Gaia      301877   Norte
3         Matosinhos      174339   Norte
4             Almada      174033  Lisboa
População de Vila Nova de Gaia: 301877
Área de Portugal: 92212

9.6 JSON INE

import requests
import os
# Ler Dados Inicial para JSON
# Indicador 0008074: Taxa de criminalida, último ano, todos os níveis geográficos, indicador 
# Categorias no SMI do Dim3: http://smi-i.ine.pt/Versao/Detalhes/902 
proxies = {
  'http': 'http://proxy.ine.pt:8080',
  'https': 'http://proxy.ine.pt:8080',
}

# os.environ['http_proxy'] = 'http://proxy.ine.pt:8080'
# os.environ['https_proxy'] = 'http://proxy.ine.pt:8080'

# Dim1=T: Dados de todos os anos
url = "https://www.ine.pt/ine/json_indicador/pindica.jsp?op=2&varcd=0008074&Dim1=T&Dim3=3&lang=PT"
print(url)
# Make an HTTP GET request to fetch the JSON data from the URL
response = requests.get(url)#, proxies=proxies)

# Verificar Resposta
# Respostas Possiveis
if response.status_code == 200:
    #Obter JSON response
    json_data = response.json()
else:
    print("Failed to fetch data. Status code:", response.status_code)
    exit()


# Verificar Tipo de Dados Devolvido e mostrar informação
if isinstance(json_data, list):
    print("JSON object is a list.")
    if json_data:
        print("Numero Registos:", len(json_data))
        #print("Registo Exemplo:", json_data[0])
        print("Tipo 1º elementos:", type(json_data[0]))
elif isinstance(json_data, dict):
    print("JSON object is a dictionary.")
    print("Keys do Object:", list(json_data.keys()))
else:
    print("Unknown JSON object type.")

    
# Obter tipo de keys no dictionary
print("Keys existentes:", list(json_data[0].keys()))

# Key com os dados

https://www.ine.pt/ine/json_indicador/pindica.jsp?op=2&varcd=0008074&Dim1=T&Dim3=3&lang=PT

JSON object is a list.
Numero Registos: 1
Tipo 1º elementos: <class 'dict'>
Keys existentes: ['IndicadorCod', 'IndicadorDsg', 'MetaInfUrl', 'DataExtracao', 'DataUltimoAtualizacao', 'UltimoPref', 'Dados', 'Sucesso']

analisar o objecto JSON devolvido

# Para poder Importar será necessário de fazer uma análise do Objeto Devolvido

# Obter Keys no Dados
# Existe um Key para Cada Ano
print("Keys existentes nos Dados:", list(json_data[0]['Dados'].keys()) )

# Ver Tipo de conteudo 2022:
print("Tipo Objecto:", type(json_data[0]['Dados']['2022']) )

# Tipo é Listagem de Dictionary's
# Ver conteudo e informação 1º elemento
print("Tipo primeiro elemento:", type(json_data[0]['Dados']['2022'][0]), 'Numero Elementos:', len(json_data[0]['Dados']['2022']) )
# Atributos de cada dictionary
print("Keys existentes no Ano:", list(json_data[0]['Dados']['2022'][0].keys()) )

Keys existentes nos Dados: ['2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019', '2020', '2021', '2022']
Tipo Objecto: <class 'list'>
Tipo primeiro elemento: <class 'dict'> Numero Elementos: 344
Keys existentes no Ano: ['geocod', 'geodsg', 'dim_3', 'dim_3_t', 'ind_string', 'valor']

mostrar ano

# Fazer um Loop por todos os anos
for ky in list(json_data[0]['Dados'].keys()):
    print(ky)

2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022

criar dataframe dos dados

import geopandas as gpd

dadosmn = r"data\geo\GPK_CAOP_MN.gpkg"
# Ler os dados do GeoPackage para um GeoDataFrame
gdfmn = gpd.read_file(dadosmn, encoding='utf-8')
print(gdfmn.info())
print(gdfmn.head())

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 308 entries, 0 to 307
Data columns (total 12 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    308 non-null    int64   
 1   DTMN        308 non-null    object  
 2   DTMNDSG     308 non-null    object  
 3   ILHA        30 non-null     object  
 4   NUTS3_02    308 non-null    object  
 5   NUTS3_15    308 non-null    object  
 6   NUTS3_24    308 non-null    object  
 7   OBJECTID_1  278 non-null    float64 
 8   NUTS1_15    278 non-null    object  
 9   NUTS2_15    278 non-null    object  
 10  DISTRITO    278 non-null    object  
 11  geometry    308 non-null    geometry
dtypes: float64(1), geometry(1), int64(1), object(9)
memory usage: 29.0+ KB
None
   OBJECTID  DTMN                DTMNDSG                ILHA NUTS3_02  \
0         1  4901                  Corvo       Ilha do Corvo      200   
1         2  4801       Lajes das Flores     Ilha das Flores      200   
2         3  4802  Santa Cruz das Flores     Ilha das Flores      200   
3         1  4201                  Lagoa  Ilha de São Miguel      200   
4         2  4202               Nordeste  Ilha de São Miguel      200   

  NUTS3_15 NUTS3_24  OBJECTID_1 NUTS1_15 NUTS2_15 DISTRITO  \
0    PT200      200         NaN     None     None     None   
1    PT200      200         NaN     None     None     None   
2    PT200      200         NaN     None     None     None   
3      200      200         NaN     None     None     None   
4      200      200         NaN     None     None     None   

                                            geometry  
0  MULTIPOLYGON (((-3463610.566 4817991.232, -346...  
1  MULTIPOLYGON (((-3479481.263 4792025.893, -347...  
2  MULTIPOLYGON (((-3474231.695 4797050.671, -347...  
3  MULTIPOLYGON (((-2845106.248 4547845.673, -284...  
4  MULTIPOLYGON (((-2812174.924 4560067.448, -281...  

import pandas as pd
# Os dados para importar dzem respeito a uma listagem de dictionary's. 

# Os dados podem ser importados a partir destas listagem com função pd.DataFrame()
# Para assegurar o tipo de dados deveria ser especificado o tipo de atributos das colunas 
columns = ["geocod", "geodsg", "dim_3", "dim_3_t", "valor"]
data_types = {"geocod": str, "geodsg": str, "dim_3": str, "dim_3_t": str, "valor": float}

# Convert the list of dictionaries to a Pandas DataFrame
df_ine = pd.DataFrame(json_data[0]['Dados']['2022'], columns=columns).astype(data_types)

# Mostrar o Resultado:
print(df_ine.info())
print(df_ine.head(8))
print('Numero registos:',len(df_ine))

# Novo Cell - filter NUTS3 (length geocod == 3):
# Alternativa - seria criar uma listagem unica a 34]
df_nuts3 = df_ine[df_ine['geocod'].str.len() == 3]
print(df_nuts3.info())

#df_nuts3['codmn'] = df_nuts3['geocod'].str[-4:]
# df_nuts3['geocod'].str[-4:]
#print(df_nuts3.head(10))

# df_mn = df_ine[df_ine['geocod'].str.len() == 7]
# print(df_mn.info())
# df_mn['codmn'] = df_mn['geocod'].str[-4:]
# print(df_mn.head(10))

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 344 entries, 0 to 343
Data columns (total 5 columns):
 #   Column   Non-Null Count  Dtype  
---  ------   --------------  -----  
 0   geocod   344 non-null    object 
 1   geodsg   344 non-null    object 
 2   dim_3    344 non-null    object 
 3   dim_3_t  344 non-null    object 
 4   valor    310 non-null    float64
dtypes: float64(1), object(4)
memory usage: 13.6+ KB
None
    geocod        geodsg dim_3                                   dim_3_t  \
0  11E0412       Vinhais     3  Furto de veículo e em veículo motorizado   
1       15       Algarve     3  Furto de veículo e em veículo motorizado   
2      150       Algarve     3  Furto de veículo e em veículo motorizado   
3  1500801     Albufeira     3  Furto de veículo e em veículo motorizado   
4  1500802      Alcoutim     3  Furto de veículo e em veículo motorizado   
5  1500803       Aljezur     3  Furto de veículo e em veículo motorizado   
6  1500804  Castro Marim     3  Furto de veículo e em veículo motorizado   
7  1500805          Faro     3  Furto de veículo e em veículo motorizado   

   valor  
0    0.0  
1    4.1  
2    4.1  
3    5.1  
4    1.6  
5   10.4  
6    3.7  
7    3.3  
Numero registos: 344
<class 'pandas.core.frame.DataFrame'>
Index: 25 entries, 2 to 342
Data columns (total 5 columns):
 #   Column   Non-Null Count  Dtype  
---  ------   --------------  -----  
 0   geocod   25 non-null     object 
 1   geodsg   25 non-null     object 
 2   dim_3    25 non-null     object 
 3   dim_3_t  25 non-null     object 
 4   valor    22 non-null     float64
dtypes: float64(1), object(4)
memory usage: 1.2+ KB
None

visualizar os dados como um mapa

# Mostrar valores unicos chave
import numpy as np
print(np.sort(df_nuts3.geocod.unique()))
print(df_nuts3[['geocod','geodsg','valor']])

# Corrigir Valores NaN
df_nuts3 = df_nuts3.fillna(0)

['111' '112' '119' '11A' '11B' '11C' '11D' '11E' '150' '16B' '16D' '16E'
 '16F' '16G' '16H' '16I' '16J' '170' '181' '184' '185' '186' '187' '200'
 '300']
    geocod                        geodsg  valor
2      150                       Algarve    4.1
22     16B                         Oeste    2.4
34     16D              Região de Aveiro    2.3
46     16E             Região de Coimbra    1.8
51     111                    Alto Minho    1.7
62     112                        Cávado    2.7
69     119                           Ave    1.6
78     11A   Área Metropolitana do Porto    4.7
79     11C                Tâmega e Sousa    1.8
87     11D                         Douro    1.3
107    16G              Viseu Dão Lafões    1.1
122    16H                   Beira Baixa    NaN
129    16I                    Médio Tejo    NaN
151    181              Alentejo Litoral    2.5
157    184                Baixo Alentejo    1.5
177    16F              Região de Leiria    2.0
196    185               Lezíria do Tejo    2.6
209    11B                   Alto Tâmega    1.1
228    186                 Alto Alentejo    1.2
240    300    Região Autónoma da Madeira    1.9
257    187              Alentejo Central    1.4
274    200    Região Autónoma dos Açores    2.3
299    16J     Beiras e Serra da Estrela    NaN
316    170  Área Metropolitana de Lisboa    3.3
342    11E      Terras de Trás-os-Montes    0.7

# Import packages
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

# Definir Figura e Axis
f, ax = plt.subplots(1, figsize=(9, 9))


# Ler os dados NUTS3
gpk = r"data\geo\GPK_NUTS3.gpkg"

# Ler os dados do GeoPackage para um GeoDataFrame
gdfnuts3 = gpd.read_file(gpk)
print(gdfnuts3.info())
# Selecionar Dados Portugal Continental:
# Fazer Seleção da NUTS1, Atributo NUTS3
# Sem seleção a area da visualização é muito grande
gdf_nuts3_sel = gdfnuts3[gdfnuts3['NUTS3'].str.startswith('1')]

print(gdf_nuts3_sel.info())

# Fazer Merge dos dados
# Fazer o Join, especificar: DF
gdf_nuts3_2 = gdf_nuts3_sel.merge(df_nuts3, left_on='NUTS3', right_on='geocod', how='left')

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 2}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf_nuts3_2.plot(column=gdf_nuts3_2.valor, 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=7, 
                      cmap='YlGn', 
                      legend=True,
                      edgecolor = 'dimgray',
                      legend_kwds  = lgnd_kwds,
                      ax=ax)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title(json_data[0]['IndicadorDsg']) # usar a designação do indicador no titulo do mapa
plt.show()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 25 entries, 0 to 24
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype   
---  ------     --------------  -----   
 0   NUTS3      25 non-null     object  
 1   NUTS3_DSG  25 non-null     object  
 2   geometry   25 non-null     geometry
dtypes: geometry(1), object(2)
memory usage: 732.0+ bytes
None
<class 'geopandas.geodataframe.GeoDataFrame'>
Index: 23 entries, 0 to 22
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype   
---  ------     --------------  -----   
 0   NUTS3      23 non-null     object  
 1   NUTS3_DSG  23 non-null     object  
 2   geometry   23 non-null     geometry
dtypes: geometry(1), object(2)
memory usage: 736.0+ bytes
None

gdfnuts3.info()
print(gdf_nuts3_sel.NUTS3.unique())
print(df_nuts3.geocod.unique())
# Valores unicos NUTS3 

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 25 entries, 0 to 24
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype   
---  ------     --------------  -----   
 0   NUTS3      25 non-null     object  
 1   NUTS3_DSG  25 non-null     object  
 2   geometry   25 non-null     geometry
dtypes: geometry(1), object(2)
memory usage: 732.0+ bytes
['111' '112' '119' '11A' '11B' '11C' '11D' '11E' '150' '16B' '16D' '16E'
 '16F' '16G' '16H' '16I' '16J' '170' '181' '184' '185' '186' '187']
['150' '16B' '16D' '16E' '111' '112' '119' '11A' '11C' '11D' '16G' '16H'
 '16I' '181' '184' '16F' '185' '11B' '186' '300' '187' '200' '16J' '170'
 '11E']

9.6.1 Exercício

Observações: - Faz a Adaptação do código seguinte para importar outro indicador ao nível de municipio ou NUTS3 - Ler Dados Indicador XXXX de um ano a escolha - Ver pagina SMI Indicadores: https://smi.ine.pt/Indicador?clear=True - Testar url antes de incluir no script - Importar para Pandas Dataframe as áreas NUTS3 - Criar mapa dos resultados

Informação URL - varcd: código de difusão - Dim1 (periódo de referência): - Ano (de acordo com portal, por exemplo S7A2016<) - Sem valores, devolvido dados último ano - T: Dados de todos os anos - Dim2: - Sem dados: retornados dados todas as geografias - Geografias separados por vírgula - Dim3=T: Informação disponível no SMI

import requests
import os
import pandas as pd

# Ler Dados Inicial para JSON
# Indicador 0008074: Taxa de criminalida, último ano, todos os níveis geográficos, indicador 
# Categorias no SMI do Dim3: http://smi-i.ine.pt/Versao/Detalhes/902 
proxies = {
  'http': 'http://proxy.ine.pt:8080',
  'https': 'http://proxy.ine.pt:8080',
}

# os.environ['http_proxy'] = 'http://proxy.ine.pt:8080'
# os.environ['https_proxy'] = 'http://proxy.ine.pt:8080'

# Dim1: Ultimo ano, Dim2: Todas as geografias 
url = r'https://www.ine.pt/ine/json_indicador/pindica.jsp?op=2&varcd=0008265&lang=PT'

# Make an HTTP GET request to fetch the JSON data from the URL
response = requests.get(url)#, proxies=proxies)

# Verificar Resposta
# Respostas Possiveis
if response.status_code == 200:
    #Obter JSON response
    json_data = response.json()
else:
    print("Failed to fetch data. Status code:", response.status_code)
    exit()

print (type(json_data))


# Verificar Tipo de Dados Devolvido e mostrar informação
if isinstance(json_data, list):
    print("JSON object is a list.")
    if json_data:
        print("Numero Registos:", len(json_data))
        #print("Registo Exemplo:", json_data[0])
        print("Tipo 1º elementos:", type(json_data[0]))
elif isinstance(json_data, dict):
    print("JSON object is a dictionary.")
    print("Keys do Object:", list(json_data.keys()))
else:
    print("Unknown JSON object type.")

    
# Obter tipo de keys no dictionary
print("Keys existentes:", list(json_data[0].keys()))

<class 'list'>
JSON object is a list.
Numero Registos: 1
Tipo 1º elementos: <class 'dict'>
Keys existentes: ['IndicadorCod', 'IndicadorDsg', 'MetaInfUrl', 'DataExtracao', 'DataUltimoAtualizacao', 'UltimoPref', 'Dados', 'Sucesso']

# Fazer um Loop por todos os anos
for ky in list(json_data[0]['Dados'].keys()):
    print(ky)

2022

# Obter Keys no Dados
# Existe um Key para Cada Ano
# Resultado json_data
print("Keys existentes nos Dados:", list(json_data[0]['Dados'].keys()) )

# Ver Tipo de conteudo 2022:
print("Tipo Objecto:", type(json_data[0]['Dados']['2022']) )

# Tipo é Listagem de Dictionary's
# Ver conteudo e informação 1º elemento
print("Tipo primeiro elemento:", type(json_data[0]['Dados']['2022'][0]), 'Numero Elementos:', len(json_data[0]['Dados']['2022']) )
# Atributos de cada dictionary
# json_data[0] = Conteudo de resposta
# json_data[0]['Dados'] = Vamos buscar os proprios dados
# Obter os dados do ano 2022 - deste conteudo podmeos criar DataFrame: json_data[0]['Dados']['2022']
print("Keys existentes no Ano:", list(json_data[0]['Dados']['2022'][0].keys()) )

Keys existentes nos Dados: ['2022']
Tipo Objecto: <class 'list'>
Tipo primeiro elemento: <class 'dict'> Numero Elementos: 344
Keys existentes no Ano: ['geocod', 'geodsg', 'ind_string', 'valor']

# criar Dataframe:
# Para assegurar o tipo de dados deveria ser especificado o tipo de atributos das colunas 
columns = ["geocod", "geodsg", "valor"]
data_types = {"geocod": str, "geodsg": str, "valor": float}

# Convert the list of dictionaries to a Pandas DataFrame
df_ine = pd.DataFrame(json_data[0]['Dados']['2022'], columns=columns).astype(data_types)
print (df_ine.head())

    geocod      geodsg  valor
0  2004901       Corvo    4.7
1  1190314      Vizela    7.5
2  1120303       Braga    7.8
3  11C1303  Felgueiras    8.1
4  11A1310     Paredes    8.2

# Mostrar os dados ao nivel de NUTS3:
# Filtragem no DF - Seleção Length 7
df_nuts3 = df_ine[df_ine['geocod'].str.len() == 7].copy()
df_nuts3['codmn'] = df_nuts3['geocod'].str[-4:]
print(df_nuts3.head())

    geocod      geodsg  valor codmn
0  2004901       Corvo    4.7  4901
1  1190314      Vizela    7.5  0314
2  1120303       Braga    7.8  0303
3  11C1303  Felgueiras    8.1  1303
4  11A1310     Paredes    8.2  1310

# ImportR gEOPackage
# Import packages
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
 
# Ler os dados CAOP
gpk = r"data\geo\GPK_CAOP_MN.gpkg"

# Ler os dados do GeoPackage para um GeoDataFrame
gdfnuts3 = gpd.read_file(gpk)
print(gdfnuts3.head())

   OBJECTID  DTMN                DTMNDSG                ILHA NUTS3_02  \
0         1  4901                  Corvo       Ilha do Corvo      200   
1         2  4801       Lajes das Flores     Ilha das Flores      200   
2         3  4802  Santa Cruz das Flores     Ilha das Flores      200   
3         1  4201                  Lagoa  Ilha de São Miguel      200   
4         2  4202               Nordeste  Ilha de São Miguel      200   

  NUTS3_15 NUTS3_24  OBJECTID_1 NUTS1_15 NUTS2_15 DISTRITO  \
0    PT200      200         NaN     None     None     None   
1    PT200      200         NaN     None     None     None   
2    PT200      200         NaN     None     None     None   
3      200      200         NaN     None     None     None   
4      200      200         NaN     None     None     None   

                                            geometry  
0  MULTIPOLYGON (((-3463610.566 4817991.232, -346...  
1  MULTIPOLYGON (((-3479481.263 4792025.893, -347...  
2  MULTIPOLYGON (((-3474231.695 4797050.671, -347...  
3  MULTIPOLYGON (((-2845106.248 4547845.673, -284...  
4  MULTIPOLYGON (((-2812174.924 4560067.448, -281...  

# Merge dos Dados:
gdf_nuts3_2 = gdfnuts3.merge(df_nuts3, left_on='DTMN', right_on='codmn', how='left')
print(gdf_nuts3_2.head())

   OBJECTID  DTMN                DTMNDSG                ILHA NUTS3_02  \
0         1  4901                  Corvo       Ilha do Corvo      200   
1         2  4801       Lajes das Flores     Ilha das Flores      200   
2         3  4802  Santa Cruz das Flores     Ilha das Flores      200   
3         1  4201                  Lagoa  Ilha de São Miguel      200   
4         2  4202               Nordeste  Ilha de São Miguel      200   

  NUTS3_15 NUTS3_24  OBJECTID_1 NUTS1_15 NUTS2_15 DISTRITO  \
0    PT200      200         NaN     None     None     None   
1    PT200      200         NaN     None     None     None   
2    PT200      200         NaN     None     None     None   
3      200      200         NaN     None     None     None   
4      200      200         NaN     None     None     None   

                                            geometry   geocod  \
0  MULTIPOLYGON (((-3463610.566 4817991.232, -346...  2004901   
1  MULTIPOLYGON (((-3479481.263 4792025.893, -347...  2004801   
2  MULTIPOLYGON (((-3474231.695 4797050.671, -347...  2004802   
3  MULTIPOLYGON (((-2845106.248 4547845.673, -284...  2004201   
4  MULTIPOLYGON (((-2812174.924 4560067.448, -281...  2004202   

                  geodsg  valor codmn  
0                  Corvo    4.7  4901  
1       Lajes das Flores   13.3  4801  
2  Santa Cruz das Flores   13.5  4802  
3                  Lagoa    8.8  4201  
4               Nordeste   17.1  4202  

# Definir Figura e Axis
f, ax = plt.subplots(1, figsize=(9, 9))

# Mostrar Dados 
# Definir Legenda 

lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 2}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf_nuts3_2.plot(column=gdf_nuts3_2.valor, 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=7, 
                      cmap='YlGn', 
                      legend=True,
                      edgecolor = 'dimgray',
                      legend_kwds  = lgnd_kwds,
                      ax=ax)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Taxa bruta de mortalidade (‰) por Local de residência (NUTS - 2013)')
plt.show()

9.7 Geocoding

API com serviços REST

API Keys google e BING (vão ser eliminadas após a formação)

• GeoCode Key BING: At0TxnfnmV0hqD99JAtRPIZfPfQarPox_JCIPgRERq-cY99c1HLvqryhnkMLwIK0
• GeoCode Key Google: AIzaSyC-tGOoI4QrYNS3AgRuzOOMb_51Gd0RTic

# definir as variaveis proxy

import os
os.environ['http_proxy'] = 'http://proxy.ine.pt:8080'
os.environ['https_proxy'] = 'http://proxy.ine.pt:8080'

obter longitude e latitudede uma morada com o Google Maps API

# Incluir Controlo de Resposta - Invalido API Key

import requests
import random, time
import pprint

proxies = {
  'http': 'http://proxy.ine.pt:8080',
  'https': 'http://proxy.ine.pt:8080',
}


API_KEY = "AIzaSyC-tGOoI4QrYNS3AgRuzOOMb_51Gd0RTic"

def geocode_address(address):
    url = "https://maps.googleapis.com/maps/api/geocode/json?address=" + address + "&key=" + API_KEY
    print(url)
    response = requests.get(url)#, proxies=proxies)
    # Mostrar Resposta JSON (para fim demonstrativos)
    print(response)
    # Atenção - ao utilizar chave errado - Response é differente     
    if response.status_code == 200:
        data = response.json()
        pp = pprint.PrettyPrinter(indent=4)
        pp.pprint(data)
        latitude = data["results"][0]["geometry"]["location"]["lat"]
        longitude = data["results"][0]["geometry"]["location"]["lng"]
        return latitude, longitude
    else:
        return None

latitude, longitude = geocode_address("Rua João Morais Barbosa 12, Lisboa")
print(latitude, longitude) 

# Para Assegurar de não ultrapassar o limite de 2 pedidos por segundo seria necessário acresentar codigo deste tipo:
time.sleep(random.uniform(0, 3)+0.1)

https://maps.googleapis.com/maps/api/geocode/json?address=Rua João Morais Barbosa 12, Lisboa&key=AIzaSyC-tGOoI4QrYNS3AgRuzOOMb_51Gd0RTic

<Response [200]>
{   'results': [   {   'address_components': [   {   'long_name': '12',
                                                     'short_name': '12',
                                                     'types': [   'street_number']},
                                                 {   'long_name': 'Rua João '
                                                                  'Morais '
                                                                  'Barbosa',
                                                     'short_name': 'R. João '
                                                                   'Morais '
                                                                   'Barbosa',
                                                     'types': ['route']},
                                                 {   'long_name': 'Lisboa',
                                                     'short_name': 'Lisboa',
                                                     'types': [   'locality',
                                                                  'political']},
                                                 {   'long_name': 'Carnide',
                                                     'short_name': 'Carnide',
                                                     'types': [   'administrative_area_level_3',
                                                                  'political']},
                                                 {   'long_name': 'Lisboa',
                                                     'short_name': 'Lisboa',
                                                     'types': [   'administrative_area_level_2',
                                                                  'political']},
                                                 {   'long_name': 'Lisboa',
                                                     'short_name': 'Lisboa',
                                                     'types': [   'administrative_area_level_1',
                                                                  'political']},
                                                 {   'long_name': 'Portugal',
                                                     'short_name': 'PT',
                                                     'types': [   'country',
                                                                  'political']},
                                                 {   'long_name': '1600-416',
                                                     'short_name': '1600-416',
                                                     'types': ['postal_code']}],
                       'formatted_address': 'R. João Morais Barbosa 12, '
                                            '1600-416 Lisboa, Portugal',
                       'geometry': {   'location': {   'lat': 38.7680256,
                                                       'lng': -9.1838724},
                                       'location_type': 'ROOFTOP',
                                       'viewport': {   'northeast': {   'lat': 38.7693255802915,
                                                                        'lng': -9.1825668697085},
                                                       'southwest': {   'lat': 38.76662761970851,
                                                                        'lng': -9.185264830291505}}},
                       'place_id': 'ChIJS_NdXMkyGQ0RlXVKqq__eac',
                       'plus_code': {   'compound_code': 'QR98+6F Lisbon, '
                                                         'Portugal',
                                        'global_code': '8CCGQR98+6F'},
                       'types': ['street_address']}],
    'status': 'OK'}
38.7680256 -9.1838724

9.7.1 Geopandas

import geopandas as gpd

geocode de uma morada

# GeoReference Morada simples:
import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import pandas as pd
import folium
# import os
# os.environ['http_proxy'] = 'http://proxy.ine.pt:8080'
# os.environ['https_proxy'] = 'http://proxy.ine.pt:8080'


# Chamar Função GeoCode
# Rua  Prof Luciano Mota vieira 42, Ponta Delgada
res_geo = gpd.tools.geocode("Rua Prof Luciano Mota vieira 42, Ponta Delgada",
                            provider = "nominatim", 
                            user_agent="Intro Geocode")

print(res_geo)

# # Atenção esta linha dá um erro caso não foi obtido nenhum resultado
# # Sera necessário fazer a validação de existencia de Geometria (utiliza Shapely)
# if (not res_geo['geometry'][0].is_empty):
#     print ('Visualizar')
#     res_geo.explore(marker_type = 'marker',edgecolor = 'black')
# else:
#     print('Nao foi obtido nenhum resultado')
# 
# print(len(res_geo))
# res_geo.explore(marker_type = 'marker',edgecolor = 'black')

                     geometry  \
0  POINT (-25.67936 37.74747)   

                                             address  
0  42, Rua Professor Luciano Mota Vieira, Ponta D...  

geocode de um ficheiro com moradas

# GeoReference Morada simples:
import geopandas as gpd
import pandas as pd
import os
from shapely.geometry import Point

# os.environ['http_proxy'] = 'http://proxy.ine.pt:8080'
# os.environ['https_proxy'] = 'http://proxy.ine.pt:8080'

# Ler Ficheiros ocm Moradas
inputfile = r"data\geo\Ensino_Nao_Superior_Amadora.xlsx"
pd_escolas = pd.read_excel(inputfile)

# Criar nova coluna morada
pd_escolas['morada2'] = pd_escolas['MORADA'].astype(str) + ' ' + pd_escolas['CTT_COD'].astype(str) + ' ' + pd_escolas['CTT_AUX'].astype(str) + ' ' + pd_escolas['LOCALIDADE'].astype(str)


# Fazer Seleção de apenas alguns registos:
pd_escolas = pd_escolas.head(10)

# Chamar Função GeoCode
try:
    res_geo = gpd.tools.geocode(pd_escolas.morada2)
except Exception as e:
    print(f"Aconteceu um erro a utilizar o geocode() de geopandas: {e}")    
    
    
#print(res_geo.info())
#print(pd_escolas.head())

print (f"Nº de Registos do ficheiro: {len(pd_escolas)}")

Nº de Registos do ficheiro: 10

selecionar os registos com geometria

from shapely.geometry import Point
import geopandas as gpd

# Mostrar Resultado:
# Selecionar Pontos com Geometria:
# Crie uma máscara booleana para identificar geometrias válidas (resultado Pandas Series)
# ~: Siginifca not (será seleccionado o inverso) - 
mask_valid_geometry = ~res_geo['geometry'].apply(lambda x: x.is_empty)

# Selecione os registros com geometria válida
res_geo_valid = res_geo[mask_valid_geometry]


print (f"Nº de Registos do ficheiro: {len(pd_escolas)}","\n",
      f"Nº de Registos resultado: {len(res_geo_valid)}")


# Mostrar a geografia obtida
res_geo_valid.explore(marker_type = 'marker',edgecolor = 'black')

Nº de Registos do ficheiro: 10 
 Nº de Registos resultado: 7

Make this Notebook Trusted to load map: File -> Trust Notebook

9.7.2 Geopy

geocode com nomination

from geopy.geocoders import Nominatim

# Morada para geocode
address = "Rua do comercio 42, vilar formoso"

# Inicialização do geocodificador com o serviço Nominatim
geolocator = Nominatim(user_agent="my_geocoder")

# Geocode da morada
# Opções para mostrar: addressdetails =True, limit = 4, extratags  = Tru
# 
location = geolocator.geocode(address)
print(location)

# Verificando o tipo de resultado JSON retornado
if location:
    print(f"Latitude: {location.latitude}, Longitude: {location.longitude}")
    print(f"Tipo de objeto JSON retornado: {type(location.raw)}")
    print("Exemplo de parte do JSON retornado:")
    print(location.raw)
else:
    print("Morada não encontrada ou geocodificação não foi possível.")

Rua do Comércio, Vilar Formoso, Almeida, Guarda, 6530-270, Portugal
Latitude: 40.6094181, Longitude: -6.8305004
Tipo de objeto JSON retornado: <class 'dict'>
Exemplo de parte do JSON retornado:
{'place_id': 252240418, 'licence': 'Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright', 'osm_type': 'way', 'osm_id': 72533753, 'lat': '40.6094181', 'lon': '-6.8305004', 'class': 'highway', 'type': 'residential', 'place_rank': 26, 'importance': 0.10000999999999993, 'addresstype': 'road', 'name': 'Rua do Comércio', 'display_name': 'Rua do Comércio, Vilar Formoso, Almeida, Guarda, 6530-270, Portugal', 'boundingbox': ['40.6067694', '40.6119459', '-6.8322906', '-6.8288417']}

geocode com google

from geopy.geocoders import GoogleV3

# Sua chave de API do Google Maps
api_key = 'AIzaSyC-tGOoI4QrYNS3AgRuzOOMb_51Gd0RTic'

# Endereço para geocode
address = "Av Antonio José Almeida, Lisboa"

# Inicialização do geocodificador com o serviço Google Maps usando a chave de API
geolocator = GoogleV3(api_key=api_key)


# Geocode da morada
location = geolocator.geocode(address)
print (location)

# Verificando os resultados
if location:
    print(f"Latitude: {location.latitude}, Longitude: {location.longitude}")
    print(location.raw)
else:
    print("Morada não encontrada ou geocodificação não foi possível.")

Av. de António José de Almeida, 1000 Lisboa, Portugal
Latitude: 38.7379621, Longitude: -9.1403984
{'address_components': [{'long_name': 'Avenida de António José de Almeida', 'short_name': 'Av. de António José de Almeida', 'types': ['route']}, {'long_name': 'Lisboa', 'short_name': 'Lisboa', 'types': ['locality', 'political']}, {'long_name': 'Lisboa', 'short_name': 'Lisboa', 'types': ['administrative_area_level_2', 'political']}, {'long_name': 'Lisboa', 'short_name': 'Lisboa', 'types': ['administrative_area_level_1', 'political']}, {'long_name': 'Portugal', 'short_name': 'PT', 'types': ['country', 'political']}, {'long_name': '1000', 'short_name': '1000', 'types': ['postal_code', 'postal_code_prefix']}], 'formatted_address': 'Av. de António José de Almeida, 1000 Lisboa, Portugal', 'geometry': {'bounds': {'northeast': {'lat': 38.7387774, 'lng': -9.136275999999999}, 'southwest': {'lat': 38.7370445, 'lng': -9.1438857}}, 'location': {'lat': 38.7379621, 'lng': -9.1403984}, 'location_type': 'GEOMETRIC_CENTER', 'viewport': {'northeast': {'lat': 38.73925993029149, 'lng': -9.136275999999999}, 'southwest': {'lat': 38.7365619697085, 'lng': -9.1438857}}}, 'place_id': 'ChIJa4H4YKEzGQ0RwerVqeaRcMY', 'types': ['route']}

geocode com bing

from geopy.geocoders import Bing

# Sua chave de API do Bing Maps
bing_api_key = 'At0TxnfnmV0hqD99JAtRPIZfPfQarPox_JCIPgRERq-cY99c1HLvqryhnkMLwIK0'

# Endereço para geocode
address = "Rua da urbanização do tanque 8, Funchal"

# Inicialização do geocodificador com o serviço Bing Maps usando a chave de API
geolocator = Bing(api_key=bing_api_key)

# Geocode da morada
location = geolocator.geocode(address)

print(location.raw,'\n')

# Verificando os resultados
if location:
    print(f"Latitude: {location.latitude}, Longitude: {location.longitude}")
else:
    print("Morada não encontrada ou geocodificação não foi possível.")

{'__type': 'Location:http://schemas.microsoft.com/search/local/ws/rest/v1', 'bbox': [32.66967721496134, -16.91852369231296, 32.67740265010269, -16.90628726787942], 'name': 'Rua da Urbanização do Tanque 8, 9050 Funchal, Portugal', 'point': {'type': 'Point', 'coordinates': [32.67353993, -16.91240548]}, 'address': {'addressLine': 'Rua da Urbanização do Tanque 8', 'adminDistrict': 'Ilha da Madeira', 'adminDistrict2': 'Funchal', 'countryRegion': 'Portugal', 'formattedAddress': 'Rua da Urbanização do Tanque 8, 9050 Funchal, Portugal', 'locality': 'Funchal', 'postalCode': '9050'}, 'confidence': 'High', 'entityType': 'Address', 'geocodePoints': [{'type': 'Point', 'coordinates': [32.67353993, -16.91240548], 'calculationMethod': 'InterpolationOffset', 'usageTypes': ['Display']}, {'type': 'Point', 'coordinates': [32.6735004, -16.91243093], 'calculationMethod': 'Interpolation', 'usageTypes': ['Route']}], 'matchCodes': ['Good']} 

Latitude: 32.67353993, Longitude: -16.91240548

importar um ficheiro inteiro (com bing)

import pandas as pd
from geopy.geocoders import Bing
import geopandas as gpd
from shapely.geometry import Point

# Seus dados
inputfile = r"data\geo\Ensino_Nao_Superior_Amadora.xlsx"
bing_api_key = 'At0TxnfnmV0hqD99JAtRPIZfPfQarPox_JCIPgRERq-cY99c1HLvqryhnkMLwIK0'

# Colunas CTT_COD e CTT_AUX são importados como numero
columns_para_string = ['CTT_COD', 'CTT_AUX']

# Leitura do arquivo Excel especificando os tipos de dados das colunas
df = pd.read_excel(inputfile, dtype={col: str for col in columns_para_string})

# Concatenando os atributos desejados para formar o endereço
df['endereco'] = df['MORADA'] + ', ' + df['CTT_COD'] + ' ' + df['CTT_AUX'] + ', ' + df['LOCALIDADE']

# Importar apenas alguns registos
df = df.head(20)

# Inicializando o geocodificador com o serviço Bing
geolocator = Bing(api_key=bing_api_key)

# Função para obter a localização e a qualidade da resposta
def get_location_info(address):
    try:
        location = geolocator.geocode(address)
        return location, location.raw['confidence']
    except:
        return None, None

# Aplicando a função para obter a localização e a qualidade da resposta
df['location_info'] = df['endereco'].apply(get_location_info)


# Extraindo as coordenadas e a qualidade da resposta para colunas separadas
# Atenção a ordem longitude (x) e latitude (y)!
df['coordinates'] = df['location_info'].apply(lambda loc: (loc[0].longitude, loc[0].latitude) if loc[0] else None)
df['quality'] = df['location_info'].apply(lambda loc: loc[1] if loc[1] else None)

# Criando o GeoDataFrame com base nas coordenadas obtidas
geometry = [Point(xy) if xy else None for xy in df['coordinates']]
# Criar gdf de resultado - com indicação do CRS
gdfBing = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")



# Corrigir para registos onde não existe GeoMetry
# Neste caso a coluna geometry is Null
mask_valid_geometry = gdfBing['geometry'].notnull()

# Selecione os registros com geometria válida
gdfBing = gdfBing[mask_valid_geometry]

print (f"Nº de Registos do ficheiro: {len(df)}","\n",
      f"Nº de Registos resultado: {len(gdfBing)}")


# Mostrar GeoDataFrame resultante
print(gdfBing.head())

print()
gdfBing.info()

Nº de Registos do ficheiro: 20 
 Nº de Registos resultado: 20
   DTCCFR                                          DSG_ESTAB  \
0  111503                JARDIM DE INFÂNCIA DA FUNDAÇÃO AFID   
1  111503  ESCOLA BÁSICA DO ALTO DO MOINHO, ZAMBUJAL, AMA...   
2  111503  ESCOLA LUIS MADUREIRA (STª CASA MISERICÓRDIA D...   
3  111501  CENTRO SOCIAL E PAROQUIAL IMACULADO CORAÇÃO DE...   
4  111503  ESCOLA BÁSICA ALMEIDA GARRETT, ALFRAGIDE, AMADORA   

                                      MORADA CTT_COD CTT_AUX LOCALIDADE  \
0         RUA QUINTA DO PARAÍSO, ALTO MOINHO    2720     502    AMADORA   
1  ESTRADA DO ZAMBUJAL BAIRRO ALTO DO MOINHO    2700     000     BURACA   
2     ESTRADA DA PORTELA - QUINTA DAS TORRES    2720     461     BURACA   
3                LARGO PADRE ADRIANO PEDRALI    2610     129  ALFRAGIDE   
4               LARGO ROTARY CLUB DA AMADORA    2720     461    AMADORA   

    LATITUDE  LONGITUDE                                           endereco  \
0  38.729557  -9.212133  RUA QUINTA DO PARAÍSO, ALTO MOINHO, 2720 502, ...   
1  38.731007  -9.215389  ESTRADA DO ZAMBUJAL BAIRRO ALTO DO MOINHO, 270...   
2  38.732370  -9.212916  ESTRADA DA PORTELA - QUINTA DAS TORRES, 2720 4...   
3  38.733139  -9.220367   LARGO PADRE ADRIANO PEDRALI, 2610 129, ALFRAGIDE   
4  38.733330  -9.214101    LARGO ROTARY CLUB DA AMADORA, 2720 461, AMADORA   

                                       location_info  \
0  ((Rua Quinta do Paraíso, Amadora, Lisboa 2610,...   
1  ((Lisbon, Portugal, (38.71221924, -9.14500046)...   
2  ((Buraca, Lisbon, Portugal, (38.74211121, -9.2...   
3  ((Largo Padre Adriano Pedrali, Alfragide, Lisb...   
4  ((Largo Rotary Club da Amadora, Amadora, Lisbo...   

                  coordinates quality                   geometry  
0  (-9.21111378, 38.73031738)    High  POINT (-9.21111 38.73032)  
1  (-9.14500046, 38.71221924)     Low  POINT (-9.14500 38.71222)  
2  (-9.20955467, 38.74211121)  Medium  POINT (-9.20955 38.74211)  
3   (-9.22027038, 38.7335836)  Medium  POINT (-9.22027 38.73358)  
4  (-9.21499638, 38.73363337)    High  POINT (-9.21500 38.73363)  

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 20 entries, 0 to 19
Data columns (total 13 columns):
 #   Column         Non-Null Count  Dtype   
---  ------         --------------  -----   
 0   DTCCFR         20 non-null     int64   
 1   DSG_ESTAB      20 non-null     object  
 2   MORADA         20 non-null     object  
 3   CTT_COD        20 non-null     object  
 4   CTT_AUX        20 non-null     object  
 5   LOCALIDADE     20 non-null     object  
 6   LATITUDE       20 non-null     float64 
 7   LONGITUDE      20 non-null     float64 
 8   endereco       20 non-null     object  
 9   location_info  20 non-null     object  
 10  coordinates    20 non-null     object  
 11  quality        20 non-null     object  
 12  geometry       20 non-null     geometry
dtypes: float64(2), geometry(1), int64(1), object(9)
memory usage: 2.2+ KB

# Exportar o resultado obtido (coluna location_info):
df['location_info'].to_csv(r'data\geo\outdfbing.txt', sep='\t')

visualizar com explore

gdfBing = gdfBing.drop(columns=['location_info'])
gdfBing.explore(marker_type = 'marker',edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

9.7.3 Dados OSMX com OSMnx

OpenStreetMap (OSM)

pesquisa de dados por lugar

import osmnx as ox
import geopandas as gpd
#import matplotlib.pyplot as plt

# Definir o lugar para qual queremos dados
place = "Porto, PT"

# Verificar Existencia Place
try:
    result = ox.geocoder.geocode(place)
    print(f"The geocoded result for {place} is: {result}")
except Exception as e:
    print(f"O place não existe: {place}")
    
tags = {"highway": "bus_stop"}

try:
    gdf_bus = ox.features_from_place(place, tags)
except Exception as e:
    print(f"Não existem elementos para este tag: {tags}")    

    
print(gdf_bus[["bench",'name','network','operator','route_ref','departures_board', 'brand' ]].head(8))

gdf_bus.explore(#column = 'cuisine',
              legend = True,
            marker_type = 'marker',
                  edgecolor = 'black')

The geocoded result for Porto, PT is: (41.1494512, -8.6107884)

                        bench                    name  network operator  \
element_type osmid                                                        
node         475347093    yes  Casa da Música (Metro)  Andante     STCP   
             492290151    yes        Pêro Vaz Caminha  Andante     STCP   
             492301138    yes        Pêro Vaz Caminha  Andante     STCP   
             609719937    yes                 Lordelo  Andante     STCP   
             987217411    NaN             Rua Firmeza  Andante     STCP   
             987217940    NaN                 Moreira  Andante     STCP   
             1403538447   yes       Hospital São João  Andante     STCP   
             1471131735   NaN             Santa Luzia  Andante     STCP   

                        route_ref departures_board brand  
element_type osmid                                        
node         475347093        NaN              NaN   NaN  
             492290151        NaN              NaN   NaN  
             492301138        NaN              NaN   NaN  
             609719937        NaN              NaN   NaN  
             987217411        NaN              NaN   NaN  
             987217940        NaN              NaN   NaN  
             1403538447       NaN              NaN   NaN  
             1471131735       NaN              NaN   NaN  

Make this Notebook Trusted to load map: File -> Trust Notebook

Pesquisa de dados por extensão de uma GDF

este exemplo faz a importação dos restaurantes existentes no OSM

import osmnx as ox
import geopandas as gpd
#import matplotlib.pyplot as plt
from shapely.geometry import box

# Por exemplo Obter Dados o Extento do GPK que utilizamos os notebooks
# Processo importar 
gpk = r'data\geo\BGRI2021_1106.gpkg'
gdf1106 = gpd.read_file(gpk,encoding='utf-8')

# Obter a extensão do GeoDataFrame
xmin, ymin, xmax, ymax = gdf1106.total_bounds

# Criar um poligono que representa a extensão
extent_polygon = box(xmin, ymin, xmax, ymax)

# Atenção O CRS dos dados do OSM é EPSG 4326 - Necesistamos de transformar o poligono para 4326

# Isto pode ser efetuado em GeoPandas (alternativa package pyproj)
# Necessário de definir o CRS par apoder fazer a projeção dos dados
extent_gdf = gpd.GeoDataFrame(geometry=[extent_polygon], crs = gdf1106.crs) 

# Mudar a projeção para CRS dos dados OSM - 4326:
extent_gdf2 = extent_gdf.to_crs('EPSG:4326')

# Esta GDF consiste de 1 registo com o poligono

# Definir Tags   
tags = {"amenity": "restaurant"}

# Necessário try e except para validar o input 
try:
    gdf_restaurants = ox.features_from_polygon(extent_gdf2['geometry'][0], tags)
except Exception as e:
    print(f"Não existem elementos para este tag: {tags}")
    

# Vai dar erro se houver problema com tags ou dados existentes    
gdf_restaurants.explore(column = 'cuisine',
        legend = True,
        marker_type = 'marker',
        edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

9.7.4 Exportar dados

gdf_restaurants.to_file(r'data\geo\osm_restaurants1106.gpkg', layer='RESTAURANTS1106', driver="GPKG")

9.7.5 Exercício

• Tentar importar mais moradas e melhorar a qualidade do endereço de input
• Importar o Ficheiro “Ensino_Nao_Superior_Amadora.xlsx” utilizando Google ou Nomantim:
  • Ver o codigo de importar utilizando Bing, com a seguinte diferença

  # Inicializar o o geocodificador com o serviço Google
  geolocator = GoogleV3(api_key=google_api_key)
  
  # Funcao get_location_info para geocodificar endereco
  def get_location_info(address):
      try:
          location = geolocator.geocode(address)
          return location, location.raw['types']
      except:
          return None, None

  • Ver o codigo de importar utilizando Bing, com a seguinte diferença

  # Inicializar o o geocodificador com o serviço Nominatim
  geolocator = Nominatim(user_agent="my_geocoder")
  
  # Funcao get_location_info para geocodificar endereco
  def get_location_info(address):
      try:
          location = geolocator.geocode(address)
          return location, location.raw['osm_type']  # Adjust according to the response structure
      except:
          return None, None

  • Ajuda Geocoders GeoPY: https://geopy.readthedocs.io/en/latest/#geocoders
• Importar as escolas de Amadora utilizando o OSMnx (tag amenity e school)
  • Ver o exemplo neste notebook
• Visualizar os diferentes Resultados obtidos:
  • É possivel de utilizar o MatplotLib para visualizar, codigo exemplo (será necessário adicionar as outras gdf

  import contextily as ctx
  from shapely.geometry import Point
  
  # Criar variáveis para a figura
  f, ax = plt.subplots(1, figsize=(9, 9))
  
  # Visualizar a GDF
  gdfBing.plot(legend = False,
                 ax = ax,
                color= 'green' )
  
  # Add basemap do contextily
  ctx.add_basemap(
      ax,
      crs=gdfGoogle.crs,
      source=ctx.providers.CartoDB.VoyagerNoLabels,
  )

Atenção: Cuidado com a quantidade de endereços a georrefenciar

9.7.5.1 Georeferenciar dados com Google

import pandas as pd
from geopy.geocoders import GoogleV3
import geopandas as gpd
from shapely.geometry import Point

# Importar os Dados
inputfile = r"data\geo\Ensino_Nao_Superior_Amadora.xlsx"
google_api_key = 'AIzaSyC-tGOoI4QrYNS3AgRuzOOMb_51Gd0RTic'

# Colunas CTT_COD e CTT_AUX são importados como números
columns_para_string = ['CTT_COD', 'CTT_AUX']

# Ler EXCEl e indicar que colunas CTT_COD e CTT_AUX são texto
df = pd.read_excel(inputfile, dtype={col: str for col in columns_para_string})

# Criar nova coluna com endereco
df['endereco'] = df['MORADA'] + ', ' + df['CTT_COD'] + ' ' + df['CTT_AUX'] + ', ' + df['LOCALIDADE']

# Importar apenas alguns registos
df = df.head(20)

# Inicializar o o geocodificador com o serviço Google
geolocator = GoogleV3(api_key=google_api_key)

# Funcao get_location_info para geocodificar endereco
def get_location_info(address):
    try:
        location = geolocator.geocode(address)
        return location, location.raw['types']
    except:
        return None, None

# Aplicando a função para obter a localização e a qualidade da resposta
df['location_info'] = df['endereco'].apply(get_location_info)

# Extraindo as coordenadas e a qualidade da resposta para colunas separadas
# Atenção a ordem longitude (x) e latitude (y)!
df['coordinates'] = df['location_info'].apply(lambda loc: (loc[0].longitude, loc[0].latitude) if loc[0] else None)
df['quality'] = df['location_info'].apply(lambda loc: loc[1] if loc[1] else None)

# Criar o GeoDataFrame com base nas coordenadas obtidas
geometry = [Point(xy) if xy else None for xy in df['coordinates']]
# Criar gdf de resultado - com indicação do CRS
gdfGoogle = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")

# Corrigir para registos onde não existe GeoMetry
# Neste caso a coluna geometry is Null
mask_valid_geometry = gdfGoogle['geometry'].notnull()

# Selecione os registros com geometria válida
gdfGoogle = gdfGoogle[mask_valid_geometry]

print (f"Nº de Registos do ficheiro: {len(df)}","\n",
      f"Nº de Registos resultado: {len(gdfGoogle)}")

# Mostrar parte do Resultado
print(gdfGoogle.head())

C:\Users\bruno.lima\AppData\Roaming\Python\Python312\site-packages\openpyxl\packaging\core.py:99: DeprecationWarning: datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).
  now = datetime.datetime.utcnow()

Nº de Registos do ficheiro: 20 
 Nº de Registos resultado: 20
   DTCCFR                                          DSG_ESTAB  \
0  111503                JARDIM DE INFÂNCIA DA FUNDAÇÃO AFID   
1  111503  ESCOLA BÁSICA DO ALTO DO MOINHO, ZAMBUJAL, AMA...   
2  111503  ESCOLA LUIS MADUREIRA (STª CASA MISERICÓRDIA D...   
3  111501  CENTRO SOCIAL E PAROQUIAL IMACULADO CORAÇÃO DE...   
4  111503  ESCOLA BÁSICA ALMEIDA GARRETT, ALFRAGIDE, AMADORA   

                                      MORADA CTT_COD CTT_AUX LOCALIDADE  \
0         RUA QUINTA DO PARAÍSO, ALTO MOINHO    2720     502    AMADORA   
1  ESTRADA DO ZAMBUJAL BAIRRO ALTO DO MOINHO    2700     000     BURACA   
2     ESTRADA DA PORTELA - QUINTA DAS TORRES    2720     461     BURACA   
3                LARGO PADRE ADRIANO PEDRALI    2610     129  ALFRAGIDE   
4               LARGO ROTARY CLUB DA AMADORA    2720     461    AMADORA   

    LATITUDE  LONGITUDE                                           endereco  \
0  38.729557  -9.212133  RUA QUINTA DO PARAÍSO, ALTO MOINHO, 2720 502, ...   
1  38.731007  -9.215389  ESTRADA DO ZAMBUJAL BAIRRO ALTO DO MOINHO, 270...   
2  38.732370  -9.212916  ESTRADA DA PORTELA - QUINTA DAS TORRES, 2720 4...   
3  38.733139  -9.220367   LARGO PADRE ADRIANO PEDRALI, 2610 129, ALFRAGIDE   
4  38.733330  -9.214101    LARGO ROTARY CLUB DA AMADORA, 2720 461, AMADORA   

                                       location_info  \
0  ((R. Q.ta do Paraíso, 2610 Amadora, Portugal, ...   
1  ((Estrada do Zambujal, 2610 Amadora, Portugal,...   
2  ((Estr. da Portela, 2610 Amadora, Portugal, (3...   
3  ((Largo Padre Adriano Pedrali, 2610-171 Amador...   
4  ((Lgo Rotary Club da Amadora, 2610-184 Amadora...   

                coordinates  quality                   geometry  
0  (-9.2120616, 38.7301915)  [route]  POINT (-9.21206 38.73019)  
1  (-9.2160293, 38.7335231)  [route]  POINT (-9.21603 38.73352)  
2  (-9.2126978, 38.7331895)  [route]  POINT (-9.21270 38.73319)  
3  (-9.2213138, 38.7332385)  [route]  POINT (-9.22131 38.73324)  
4   (-9.2154748, 38.733811)  [route]  POINT (-9.21547 38.73381)  

# Apagar atributo location_info 
gdfGoogle = gdfGoogle.drop(columns=['location_info'])
gdfGoogle.explore(marker_type = 'marker',edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

9.7.5.2 Importar dados Nominatim

import pandas as pd
from geopy.geocoders import Nominatim
import geopandas as gpd
from shapely.geometry import Point

# Importar os Dados
inputfile = r"data\geo\Ensino_Nao_Superior_Amadora.xlsx"

# Colunas CTT_COD e CTT_AUX são importados como números
columns_para_string = ['CTT_COD', 'CTT_AUX']

# Ler EXCEl e indicar que colunas CTT_COD e CTT_AUX são texto
df = pd.read_excel(inputfile, dtype={col: str for col in columns_para_string})

# Criar nova coluna com endereco
df['endereco'] = df['MORADA'] + ', ' + df['CTT_COD'] + ' ' + df['CTT_AUX'] + ', ' + df['LOCALIDADE']

# Importar apenas alguns registos
df = df.head(20)

# Inicializar o o geocodificador com o serviço Nominatim
geolocator = Nominatim(user_agent="my_geocoder")

# Funcao get_location_info para geocodificar endereco
def get_location_info(address):
    try:
        location = geolocator.geocode(address)
        return location, location.raw['osm_type']  # Adjust according to the response structure
    except:
        return None, None

# Aplicando a função para obter a localização e a qualidade da resposta
df['location_info'] = df['endereco'].apply(get_location_info)

# Extraindo as coordenadas e a qualidade da resposta para colunas separadas
# Atenção a ordem longitude (x) e latitude (y)!
df['coordinates'] = df['location_info'].apply(lambda loc: (loc[0].longitude, loc[0].latitude) if loc[0] else None)
df['quality'] = df['location_info'].apply(lambda loc: loc[1] if loc[1] else None)

# Criar o GeoDataFrame com base nas coordenadas obtidas
geometry = [Point(xy) if xy else None for xy in df['coordinates']]
# Criar gdf de resultado - com indicação do CRS
gdfNominatim = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")

# Corrigir para registos onde não existe GeoMetry
# Neste caso a coluna geometry is Null
mask_valid_geometry = gdfNominatim['geometry'].notnull()

# Selecione os registros com geometria válida
gdfNominatim = gdfNominatim[mask_valid_geometry]

print (f"Nº de Registos do ficheiro: {len(df)}","\n",
      f"Nº de Registos resultado: {len(gdfNominatim)}")


# Mostrar GDf Resultado
print(gdfNominatim.head())

C:\Users\bruno.lima\AppData\Roaming\Python\Python312\site-packages\openpyxl\packaging\core.py:99: DeprecationWarning: datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).
  now = datetime.datetime.utcnow()

Nº de Registos do ficheiro: 20 
 Nº de Registos resultado: 5
    DTCCFR                                          DSG_ESTAB  \
4   111503  ESCOLA BÁSICA ALMEIDA GARRETT, ALFRAGIDE, AMADORA   
11  111501  ESCOLA BÁSICA DA QUINTA GRANDE, ALFRAGIDE, AMA...   
13  111503        ESCOLA BÁSICA ALICE VIEIRA, BURACA, AMADORA   
14  111504  ESCOLA BÁSICA PROF. PEDRO D¿OREY DA CUNHA, DAM...   
17  111503  CENTRO INFANTIL DE SÃO GERARDO - CENTRO SOCIAL...   

                                               MORADA CTT_COD CTT_AUX  \
4                        LARGO ROTARY CLUB DA AMADORA    2720     461   
11                                AV. DAS LARANJEIRAS    2720     334   
13                            R. PROF. DR. EGAS MONIZ    2610     150   
14                              R. BERNARDINO MACHADO    2720     066   
17  RUA DE SANTA FILOMENA - ALTO DA COVA DA MOURA ...    2610     210   

   LOCALIDADE   LATITUDE  LONGITUDE  \
4     AMADORA  38.733330  -9.214101   
11  ALFRAGIDE  38.739212  -9.221827   
13     BURACA  38.742123  -9.210683   
14     DAMAIA  38.742447  -9.219668   
17    AMADORA  38.744495  -9.211735   

                                             endereco  \
4     LARGO ROTARY CLUB DA AMADORA, 2720 461, AMADORA   
11           AV. DAS LARANJEIRAS, 2720 334, ALFRAGIDE   
13          R. PROF. DR. EGAS MONIZ, 2610 150, BURACA   
14            R. BERNARDINO MACHADO, 2720 066, DAMAIA   
17  RUA DE SANTA FILOMENA - ALTO DA COVA DA MOURA ...   

                                        location_info  \
4   ((Largo Rotary Club da Amadora, Alfragide, Ama...   
11  ((Avenida das Laranjeiras, Quinta Grande, Alfr...   
13  ((Rua Professor Doutor Egas Moniz, Alto da Cov...   
14  ((Rua Bernardino Machado, Alto da Cova da Mour...   
17  ((Rua de Santa Filomena, Alto da Cova da Moura...   

                 coordinates quality                   geometry  
4   (-9.2150865, 38.7338984)     way  POINT (-9.21509 38.73390)  
11    (-9.2216721, 38.73717)     way  POINT (-9.22167 38.73717)  
13  (-9.2115064, 38.7416111)     way  POINT (-9.21151 38.74161)  
14  (-9.2189439, 38.7432611)     way  POINT (-9.21894 38.74326)  
17  (-9.2117563, 38.7443726)     way  POINT (-9.21176 38.74437)  

# Google não tem problema com o atributo location_info
gdfNominatim = gdfNominatim.drop(columns=['location_info'])
gdfNominatim.explore(marker_type = 'marker',edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

9.7.5.3 Importar dados OSM

import osmnx as ox
import geopandas as gpd
import matplotlib.pyplot as plt

# Definir o lugar para qual queremos dados
place = "Amadora, PT"

# Verificar Existencia Place
try:
    result = ox.geocoder.geocode(place)
    print(f"The geocoded result for {place} is: {result}")
except Exception as e:
    print(f"O place não existe: {place}")
    
tags = {"amenity": "school"}

try:
    gdf_school = ox.features_from_place(place, tags)
except Exception as e:
    print(f"Não existem elementos para este tag: {tags}")    

    
gdf_school.explore(marker_type = 'marker',
                  edgecolor = 'black')

The geocoded result for Amadora, PT is: (38.758959, -9.2365233)

Make this Notebook Trusted to load map: File -> Trust Notebook

mostrar todos os resultados

import contextily as ctx
from shapely.geometry import Point

# Criar variáveis para a figura
f, ax = plt.subplots(1, figsize=(9, 9))


# Visualizar a GDF
gdfBing.plot(legend = False,
               ax = ax,
              color= 'green' )

# Add basemap do contextily
ctx.add_basemap(
    ax,
    crs=gdfGoogle.crs,
    source=ctx.providers.CartoDB.VoyagerNoLabels,
)



# Visualizar a GDF
gdfNominatim.plot(legend = False,
               ax = ax,
              color= 'purple' )

# Visualizar a GDF
gdfGoogle.plot(legend = False,
               ax = ax,
              color= 'red' )


# Visualizar a GDF
gdf_school.plot(legend = False,
               ax = ax,
              color= 'blue' )

# Add basemap do contextily
ctx.add_basemap(
    ax,
    crs=gdfGoogle.crs,
    source=ctx.providers.CartoDB.VoyagerNoLabels,
)


ax.set_axis_off()